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Abstract 

Models  are  developed  for  a  state-of-the-art  time-domain  ship  motion  program  to  predict  ship 
motions  during  flooding  and  green  water  on  deck  events.  Water  mass  from  the  flooding  and 
green  water  is  incorporated  into  the  dynamic  equations  of  motion  using  time-dependent  mass 
and  moment  of  inertia  terms. 

Green  water  on  deck  includes  three  subproblems:  the  problem  of  water  shipping  on  deck, 
the  problem  of  motion  of  water  trapped  on  the  deck,  and  the  problem  of  water  escaping  off 
the  deck.  This  research  looks  at  the  first  two  suproblems,  both  of  which  involve  shallow  water 
wave  theory.  Glimms  method,  also  called  the  Random  Choice  Method,  and  the  Flux  Difference 
Splitting  Method  are  both  investigated  as  solution  techniques  for  the  motion  of  water  on  deck. 

This  work  provides  a  tool  to  estimate  ship  damaged  stability  and  examine  the  effects  of 
progressive  flooding. 

Thesis  Supervisor:  Dick  K.P.  Yue 

Title:  Professor  of  Hydrodynamics  and  Ocean  Engineering 
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Chapter  1 

Introduction 


1.1  Background 

There  are  many  hazards  to  ships  that  can  result  in  hull  damage  and  subsequent  flooding. 
Depending  on  the  extent  of  a  ship’s  damaged  condition,  flooding  may  cause  a  loss  of  buoyancy, 
a  loss  of  transverse  stability,  and  significant  changes  in  trim  and  list.  Adverse  buoyancy  and 
trim  conditions  can  lead  to  sinking  by  foundering,  while  the  loss  of  transverse  stability  can  lead 
to  capsizing.  Significant  trim  and  list  changes  may  also  result  in  water  on  the  weatherdeck  due 
to  shipping  water  as  freeboard  is  lost.  The  water  on  deck,  often  referred  to  as  green  water  or 
the  green  water  problem,  can  further  harm  a  damaged  ship’s  stability  condition  and  also  affect 
the  main  hull  girder  loads,  and  deck  and  superstructure  loads. 

The  state  of  stability,  list,  and  trim  in  a  damaged  ship  is  dynamic;  it  varies  over  time  as 
the  flooding  event  progresses  and  also  depends  strongly  on  environmental  conditions  such  as 
sea  state  and  wind.  Current  naval  standards,  however,  take  a  static  approach  in  specifying 
stability  requirements  for  a  damaged  ship.  For  example,  the  naval  standard  DDS-079,  reference 
[5],  requires  that  stability  be  analyzed  on  the  equilibrium  position  of  the  damaged  ship  based 
purely  on  static  geometry  after  the  flooding  event  is  complete.  This  analysis  is  similar  in 
many  respects  to  intact  stability  calculations  except  with  characteristics  such  as  metacenter, 
center  of  gravity,  and  righting  arm  curves  adjusted  due  to  the  weight  of  water  in  the  flooded 
compartments.  Reference  [5]  makes  use  of  wind  speed  and  wave  height  for  damaged  stability 
analysis,  but  these  environmental  conditions  are  also  applied  to  the  analysis  in  a  static  sense 
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through  applying  steady  wind  healing  moments  and  placing  limitations  on  bulkhead  opening 
locations.  Reference  [9]  refers  to  such  bulkhead  locations  as  “V-lines.” 

In  [31]  Surko  points  out  limitations  due  to  the  static  approach  of  current  damaged  stability 
analysis  procedure  and  criteria.  Of  these  limitations,  two  are  becoming  more  salient  as  the  US 
Navy  shifts  to  performance  based  requirements.  First,  in  1987  the  Chief  of  Naval  Operations 
(CNO)  [24]  endorsed  a  series  of  operational  characteristics  to  be  incorporated  into  surface 
combatants  of  the  year  2010.  Included  in  these  characteristics  is  that  a  ship  has  the  capability 
to  fight,  even  though  it  may  have  sustained  hull  damage  and  be  flooded,  with  whatever  weapons 
systems  are  available.  To  assess  whether  a  ship  could  employ  weapons  while  fighting  hurt  would 
require  analyzing  ship  motions  which  is  well  beyond  the  scope  of  reference  [5]  procedures.  In  fact 
references  [14],  [15],  and  [16]  report  there  is  no  information  in  the  literature  and  no  appropriate 
computer  prediction  tools  to  assess  ship  motion  performance  of  partially  flooded  or  flooding 
ships  in  waves  and  wind.  The  second  significant  limitation  in  current  damaged  stability  criteria 
pointed  out  by  Surko  is  that  moderate  wind  and  sea  conditions  are  assumed.  In  reference  [31] 
Surko  shows  there  is  a  considerable  probability  of  experiencing  wave  action  that  exceeds  the 
moderate  8  foot  wave  height  assumed  in  reference  [5]. 

As  a  first  step  in  addressing  these  limitations,  model  tests  have  been  performed  to  assess 
the  dynamic  stability  of  current  fleet  combatants  in  a  damaged  condition.  These  tests  were 
reported  in  references  [14],  [15],  and  [16].  The  term  ”  dynamic  stability’  in  these  model  tests 
is  meant  in  its  true  sense:  actual  ship  motions  and  ability  to  withstand  sinking  under  a  variety 
of  environmental  and  flooding  conditions. 

Model  testing  can  be  costly  due  to  production  of  scale  models  and  the  use  of  large  labora¬ 
tory  facilities  for  the  experiments.  Also,  the  time  requirements  to  prepare  and  conduct  model 
tests  make  it  difficult  to  use  testing  early  in  the  design  process  to  predict  damaged  dynamic 
stability  for  immature  designs  and  design  variants.  Development  of  damaged  stability  computer 
prediction  tools,  especially  for  early  in  the  design  process,  would  be  ideal  as  a  supplement  or 
replacement  to  model  testing. 

Computational  fluid  dynamic  (CFD)  codes  that  predict  ship  motion  would  provide  a  good 
foundation  for  development  of  damaged  stability  prediction  tools.  Of  the  two  general  categories 
of  CFD  codes  that  predict  ship  motion,  frequency  domain  and  time  domain,  the  time  domain 
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approach  is  better  suited  for  damaged  stability  analysis.  Frequency  domain  codes  only  consider 
the  mean  underwater  hull  form  and  linearize  by  assuming  small  wave  and  motion  amplitudes. 
The  linear  prediction  would  breakdown  under  high  sea  states  and  large  amplitude  responses 
that  need  to  be  considered  for  a  damaged  stability  analysis.  On  the  other  hand,  state  of  the 
art  CFD  codes  that  predict  wave-induced  ship  motions  and  loads  in  the  time  domain  solve 
the  non-linear  three-dimensional  ship  motion  problem  and  can  handle  large  wave  and  motion 
amplitudes. 

Also,  inherent  in  use  of  a  time  domain  code  as  a  damage  stability  prediction  tool  is  the 
ability  to  predict  motions  during  the  entire  flooding  event.  Such  a  tool  would  be  useful 
at  assessing  the  effects  of  progressive  flooding.  Progressive  flooding  occurs  when  water  in 
flooded  compartments  floods  into  adjacent  compartments  by  overflowing  watertight  bulkheads 
or  leaking  through  damaged  bulkheads.  The  R.M.S.  Titanic  sank  as  a  result  of  progressive 
flooding  which  flooded  compartments  beyond  those  originally  opened  to  the  sea  by  the  iceberg- 
caused  damage.  Progressive  flooding  is  of  special  concern  in  warships  where  hull  damage  from 
combat  is  likely  to  cause  the  watertight  bulkheads  surrounding  the  affected  compartments  to 
suffer  some  damage  from  shock  or  fragmentation.  The  US  Navy  has  an  interest  in  progressive 
flooding  but  the  published  work  to  date,  an  example  of  which  is  in  reference  [2],  has  been  simple 
quasistatic  models  that  do  little  more  than  determine  the  damaged  ship  hydrostatic  position 
throughout  the  progressive  flooding  event. 

1.2  Research  Objectives 

This  research  investigates  the  addition  of  a  compartment  flooding  model  and  green  water  model 
to  a  CFD  code  that  predicts  ship  motions  in  the  time  domain  so  that  it  can  be  used  as  a 
damaged  stability  prediction  tool.  There  is  no  effort  made  by  the  author  to  perform  damaged 
stability  analysis.  The  specific  CFD  code  used  for  these  purposes  is  the  Large  Amplitude 
Motions  Program  (LAMP)  developed  by  the  Ship  Technology  Division  of  Science  Application 
International  Corporation  (SAIC).  The  theory  and  some  results  of  the  LAMP  code  have  been 
presented  in  several  papers  including  references  [7],  [26],  and  [19].  A  brief  review  of  the  theory 
and  formulations  of  LAMP  is  given  in  Chapter  2. 
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Green  water  and  compartment  flooding  can  be  considered  as  events  that,  at  each  instant, 
are  part  of  the  ship  and  change  the  total  rigid  body  mass  and  mass  distribution.  Sloshing  and 
water  motion  will  also  affect  the  mass  distribution.  Both  events  are  fundamentally  the  same 
process  that  can  be  modeled  as  a  time-rate-of-change  of  ship  mass  in  the  rigid  body  motion 
problem.  The  approach  in  this  thesis,  then,  is  to  calculate  the  affect  on  ship  motions  from  green 
water  and  flooding  by  incorporating  time-dependent  mass  and  mass  moment  of  inertia  into  the 
LAMP  dynamic  equations  of  motion  solver. 

The  green  water  problem  includes  three  subproblems:  water  shipping;  motion  of  water  on 
deck;  and  water  escaping  off  deck.  This  research  looks  at  the  first  two  subproblems  in  some 
detail.  Water  escaping  is  treated  by  simply  letting  water  fall  off  the  weatherdeck  edges. 

The  water  motion  on  deck  subproblem  involves  shallow  water  wave  theory.  There  are 
several  solution  techniques  that  have  been  developed  to  solve  the  shallow  water  wave  problem. 
Two  of  the  techniques,  Glimms  method  (also  called  the  Random  Choice  Method)  and  the  Flux 
Difference  Splitting  Method,  are  robust  in  the  sense  that  they  can  handle  dicontinuities  such  as 
shocks  and  bores  in  the  solution.  This  research  investigates  implementation  of  both  solution 
techniques.  As  a  result,  the  flux  difference  splitting  method  is  selected  as  the  solution  technique 

4 

for  shallow  water  flow  in  the  green  water  model. 
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Chapter  2 


LAMP  Description  and 
Development  of  Equations  of  Motion 


2.1  LAMP  Description 

This  description  of  LAMP  is  primarily  based  on  information  from  reference  [21].  LAMP  com¬ 
putes  a  time  domain  solution  for  a  general  three-dimensional  body  floating  on  a  free  surface. 
Six  degree-of-freedom  motions  are  permitted.  LAMP  obtains  a  potential  flow  solution  to  the 
body-wave  interaction  problem  using  the  boundary-element  (or  panel)  method  where  the  sub¬ 
merged  body  surface  is  divided  into  a  number  of  panels.  The  incoming  waves  can  take  any 
form.  At  each  time  step  the  hydrodynamic  pressure  forces  on  the  hull,  which  are  computed 
from  the  complete  velocity  potential  solution,  are  combined  with  body  forces  and  any  external 
forces  to  solve  the  equations  of  motion.  The  hull  pressure  forces  may  also  be  used  to  calculate 
hull  bending  and  torsional  moments  and  shear  forces. 

In  order  to  balance  computation  requirements  with  physics  correctness  and  complexity, 
LAMP  has  three  methods  of  calculation.  The  user  selects  a  specific  calculation  method  for  a 
LAMP  run  through  control  variables  specified  in  the  input.  The  LAMP  calculation  methods 
are  compared  in  Table  2.1. 
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Method 

Hydrodynamic,  Restoring,  and  Froude-Krylov  Wave  Forces 

LAMP-1 

Free  Surface  Boundary  Condition  on  Mean  Water  Surface 

3-D  Linear  Hydrodynamics 

Linear  Hydrostatic  Restoring  and  Froude-Krylov  Wave  Forces 

LAMP-2 

Free  Surface  Boundary  Condition  on  Mean  Water  Surface 

3-D  Linear  Hydrodynamics 

Nonlinear  Hydrostatic  Restoring  and  Froude-Krylov  Wave  Forces 

LAMP-4 

Free  Surface  Boundary  Condition  on  Incident  Water  Surface 

3-D  Nonlinear  Hydrodynamics 

Nonlinear  Hydrostatic  Restoring  and  Froude-Krylov  Wave  Forces 

Table  2.1:  LAMP  Calculation  Methods  and  Description 


The  LAMP-4  method  is  the  complete  large-amplitude  method  where  the  3-D  velocity  po¬ 
tential  is  computed  with  the  linearized  free-surface  condition  satisfied  on  the  surface  of  the 
incident  wave.  Both  the  hydrodynamic  and  hydrostatic  pressure  are  computed  over  the  in¬ 
stantaneous  hull  surface  below  the  incident  wave  surface.  The  incoming  wave  slope  must  be 
small.  Small  slope  generally  indicates  that  the  wave  height  is  one  order  of  magnitude  less  than 
the  wavelength.  LAMP-4  has  large  computational  requirements  and  has  traditionally  been 
run  on  a  supercomputer.  However,  the  mixed-source  formulation  now  used  in  LAMP  to  solve 
the  potential  flow  problem  provides  enough  computational  savings  for  LAMP-4  to  be  run  on  a 
workstation.  The  mixed  source  formulation  is  discussed  later  in  this  chapter. 

The  LAMP-2  method  is  an  approximate  nonlinear  method  which  retains  many  of  the  ad¬ 
vantages  of  both  LAMP-1  and  LAMP-4.  It  uses  a  linear  3-D  approach  like  LAMP-1,  where 
the  potential  flow  problem  is  solved  over  the  mean  body  boundary  position,  to  compute  the 
hydrodynamic  (radiation,  diffraction,  and  forward  speed)  part  of  the  pressure  forces.  The 
hydrostatic  restoring  and  Froude-Krylov  wave  forces  are  calculated  on  the  portion  of  the  ship 
beneath  the  incident  wave  surface.  The  requirements  for  computer  resources  are  about  the 
same  as  LAMP-1.  Note  that  the  LAMP-2  and  LAMP-4  nonlinear  methods  are  based  on  the 
approach  that  both  the  ship  motions  and  the  waves  may  have  large  amplitudes. 

The  LAMP-1  method  is  the  linearized  version  of  the  LAMP-4  method  with  the  free  surface 
boundary  conditions  satisfied  on  the  undisturbed  free  surface  location.  The  linear  hydrostatic 
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restoring  forces  are  computed,  from  waterplane  quantities  while  the  Froude-Krylov  wave  forces 
are  calculated  with  pressure  below  the  undisturbed  free  surface  location.  Like  LAMP-2,  the 
mean  body  boundary  position  is  used  for  the  potential  flow  problem. 

LAMP  uses  two  approaches  toward  solving  the  hydrodynamic  problem  for  the  potential 
function  $(i)  at  each  time  step:  a  direct  solution  of  the  hydrodynamic  potentials  in  the  time 
domain  and  a  solution  using  pre-computed  impulse  response  functions.  Both  solutions  are  based 
upon  a  mixed-source  formulation  that  is  briefly  described  in  the  remainder  of  this  section. 

In  the  mixed-source  formulation,  both  the  Rankine  source  and  the  transient  Green  function 
are  used.  The  fluid  domain  is  divided  into  an  inner  domain  I  and  an  outer  domain  II  as  shown 
in  figure  2-l.The  inner  domain  is  enclosed  by  the  wetted  body  surface  Sb,  a  local  portion  of  the 
free  surface  Sf,  and  the  matching  surface  Sm.  The  free  surface  Sf  intersects  the  body  surface 
and  is  truncated  by  the  matching  surface  iS^-at  the  water  fine  IVn,  ■  The  outer  domain  is  the 
rest  of  the  fluid  region  enclosed  by  Sm,  an  imaginary  surface  S^,  and  the  remaining  free  surface 
intersected  by  S'TO.and  Soo- 


Figure  2-1:  Domain  Definitions  in  the  LAMP  Mixed-Source  Formulation 
The  fluid  motion  is  described  by  a  velocity  potential, 

$T(~x,t)  =  +  <&('•?,  t)  (2.1) 
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where  is  the  incident  wave  potential  and  $  is  the  total  disturbance  potential  due  to  the 
presence  of  the  ship,  x  is  a  position  vector  and  t  is  time.  In  the  inner  domain  7,  the  initial 
boundary  value  problem  for  $  =  $>/  can  be  expressed  as, 

V24>j  =  0  in  7  (2.2) 


The  inner  domain  potential  must  satisfy  the  free  surface  and  body  boundary  conditions.  The 
free  surface  boundary  condition  is  linearized  in  all  three  formulations,  such  that 


**L  +  g°*L  =  0 

m?  9  dz 


sf(t),t>  o 


where  g  is  the  gravitational  acceleration.  The  body  boundary  condition  is  next  applied  on  the 
instantaneous  underwater  body  for  LAMP-4  and  the  mean  underwater  body  for  LAMP-1  and 
LAMP-2, 

„„  Sb(t),t  >  0  (2.4) 

d n  on 

- ^  , 

where  it  is  a  unit  normal  vector  to  the  body  out  of  the  fluid  and  V  n  is  the  instantaneous 
body  velocity  in  the  normal  direction.  S^(t)  is  constant  for  LAMP-1  and  LAMP-2.  Finally, 
the  initial  conditions  require  a  zero  disturbance  potential  on  the  free  surface  at  t  0, 


*,=  ®2i  =  o 
1  at 


The  corresponding  boundary  integral  equation  in  terms  of  the  Rankine  source  is, 

2tt$/(P)+  [  ($IGn-®inG)dS  =  0  (2.6) 

J  Sj 

where  G  =  1/r  =  1/  \P  -  Q\.  P  =  (x,y,z)  and  Q  =  (£,??,  C)  are  the  field  point  and  source 
point  on  Sj  =  Sf  U  Si,  U  Sm. 

In  the  outer  domain  II,  the  initial-value  boundary  problem  for  $  =  $//  can  be  written  as, 


V2$//  =  0  in  II 
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V$jj  — *•  0  at  co 

$//  =  ^  =  0  at  *  =  0  (2.9) 

ot 

The  corresponding  boundary  integral  equation  in  terms  of  the  Rankine  source  is, 

27 K*n(P)  +  [  (<S>nG°n  -  <t>unG°)dS  =  M(P,  t )  (2.10) 

J  Sm 

where  the  memory  fimction  M(P,  t)  is  defined  as, 

M(P, t)=  f  dr  (  [  ($nGfTn  -  <S>IInGfT)dS  +  -  [  ( <Z>iiGfTT  -  $//r GfT)VNdL\  (2.11) 

Jo  kJSm  9  > 

where  rm  is  the  water  line  of  the  matching  surface,  V/v  is  the  normal  velocity  of  rm,  and  G° 
and  Gf  are  associated  with  the  transient  Green  function.  Reference  [29]  provides  a  detailed 
description  of  the  transient  Green  function. 

The  matching  surface  is  treated  as  a  control  surface  and  moves  with  the  body.  To  complete 
the  problem  statement,  the  matching  conditions  require  that  the  total  disturbance  velocity 
potential  and  the  normal  velocity  across  the  matching  surface  are  continuous,  thereby  producing 


=  $//  on  Sm 


(2.12) 


5$/  _  d<&u 
dn  dn 


(2.13) 


The  solution  is  obtained  at  each  time  step.  Using  the  panel  method,  the  above  equations 
are  used  to  solve  for  *1' /  on  S&,  on  Sj ,  and  ‘l’/  and  ^n[  on  Sm.  Bernoulli’s  equation  is  used 

to  compute  the  pressure  on  the  hull  surface,  which  is  integrated  to  get  the  hydrodynamic  forces 
on  the  ship.  Then  the  linearized  free  surface  boundary  condition  can  be  used  in  domain  I  to 
integrate  in  time  and  update  the  values  of  the  total  disturbance  wave  elevation  and  <&/  at  the 


next  time  step. 


16 


2.2  LAMP  Rigid  Body  Dynamics 

This  section  outlines  the  solution  in  LAMP  to  the  rigid  body  dynamic  problem.  Several 
coordinate  systems,  illustrated  in  figure  2-2,  are  used  to  describe  the  six-degree-of-freedom 
motion  of  a  ship  in  a  seaway.  The  global  system  ,Og,  is  fixed  on  earth.  A  second  system  is 
the  local  system,  Of,  which  is  fixed  at  the  ship’s  center  of  mass,  eg,  and  rotates  with  the  ship. 
The  relation  between  these  two  systems  is  by  the  position  vector,  R,  and  a  set  of  euler  angles, 
=  (<Ev,©,  vp),  measured  in  the  global  system  and  following  the  sequence  of  rotation  S',  0, 
and  4>r,  respectively,  fromOg  to  O/.  The  angles  f2  can  be  thought  of  in  common  terms:  ^  is 


Figure  2-2:  Coordinate  Systems  for  LAMP  Dynamic  Solver 

the  ship’s  yaw  ,  ©  is  the  ship’s  pitch,  and  3>r  is  the  ship’s  roll.  The  matrix  L  is  the  euler  angle 
transformation  matrix  between  Og  and.O/. 


L  = 


cos  'F  cos  ©  sin  'I'  cos  ©  —  sin  © 

—  sin  ^  cos  $r  +  cos  $  sin  ©  sin  3>r  cos  cos  <f>r  +  sin  'I'  sin  0  sin  <f>r  cos©  sin  $r 
sin^sin^r  +  cos®  sin©  cos  $r  —  cos  ®  sin  3V  +  sin  ®  sin  ©  cos  <E>r  cos©  cos  $r 


(2.14) 

A  third  coordinate  system,  Og',has  the  same  orientation  as  Og  but  is  initially  centered  at  eg 
and  moves  with  steady  ship  speed,  Uship.  There  are  other  coordinate  systems  used  in  LAMP 
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to  define  the  rigid  body  geometry,  for  example  the  input  and  initial  static  systems,  but  these 
are  not  necessary  in  describing  the  solution  to  the  dynamic  equations  of  motion. 

Velocities  for  the  dynamic  problem  are  defined  as  follows.  All  linear  velocities  are  referred 
to  the  global  coordinate  system  unless  indicated.  Vo  is  the  velocity  of  a  point  on  the  rigid 
body,  or  extended  rigid  body,  that  coincides  with  Og  at  time  f;  V  is  the  rigid  body  velocity  at 
eg;  Vg  is  the  rigid  body  velocity  at  eg  referred  to  the  0'g  system,  c o  is  the  absolute  angular 
velocity  and  .  Tot  is  the  angular  velocity  in  local  system,  O/. 

The  velocities  are  related  by 


V=Vo  +  #xR 

Vg  =  Vo  —  U  ship 


Tj/  =  L~uj 

rrt d  .  — 'J1  d  y 

u t  —  L  ojf  and  —a ;  =  L  —a;/ 

dt  at 


The  rate  of  change  of  It  in  terms  of  the  angular  velocity  of  the  ship,  Tot,  is  given  by 


(2.15) 

(2.16) 

(2.17) 

(2.18) 
(2.19) 


(2.20) 


where  [W]  is  defined  as 


1  sin  sin©/ 
0  cos  <Ev 


iin©/cos©  cos  $r  sin 0/ cos  © 
os<E>r  —  sin  <hr 


0  sin <Jv/ cos©  cos  $r/ cos  © 


(2.21) 


To  determine  R  and  it  the  dynamic  equations  for  the  rigid  body  motion  must  be  solved. 
The  equation  of  motion  for  translation  can  be  written  as 


fb  =  —  ((Mo)(V„  -  Uship )) 


(2.22) 


where  fb  is  the  total  force  acting  on  the  rigid  body  at  eg  and  Mo  is  the  rigid  body  mass  matrix. 
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The  equation  of  motion  for  rotation  about  eg  in  the  local  coordinate  system  can  be  written  as 

Mo  =  LT I/j-(  w7)  +LT {~u'  x  VStl)  (2.23) 

where  Mo  is  the  total  moment  acting  on  the  rigid  body  about  eg  in  the  Og  or  Og'  coordinate 
system  and  1/  is  the  rigid  body  mass  moment  of  inertia  tensor  in  the  local  coordinate  system. 
Fo  and  Mo  are  calculated  from  the  instantaneous  total  force  including  hydrodynamics  and 
external  contributions.  Equations  2.22  and  2.23  can  be  formulated  into  a  coupled  system  of 
equations 

r  1  J  I  v*  I  r  1 

(2.24) 


1 

d 

Vo 

1 

E 

dt 

~ul 

9 

Combining  equations  2.24,  2.16,  and  2.20  yields 


d_ 

dt 


V0 

p 

_ V 

a?/ 

R 

_ y 

ft 

- 

E 


-l  r 


V 

[w\i3i 


(2.25) 


Equation  2.25  is  solved  by  the  fourth  order  Runge-Kutta  method. 


2.3  LAMP  Rigid  Body  Dynamics  With  Time-Dependent  Mass 

The  solution  to  the  rigid  body  dynamic  problem  requires  some  modification  to  account  for 
time-dependent  mass  and  mass  moment  of  inertia  tensor  due  to  water  added  to  a  ship  from 
flooding  or  shipping  water.  Also,  water  motion  will  change  the.  mass  moment  of  inertia  tensor 
over  time.  This  section  provides  a  detailed  formulation  of  the  rigid  body  dynamics  solution 
with  time-dependent  mass  and  mass  moment  of  inertia  tensor. 

2.3.1  Infinite  Frequency  Added  Mass  and  Moment  of  Inertia 

The  solution  to  the  rigid  body  dynamic  equations  of  motion  in  LAMP  makes  use  of  the  ship’s 
infinite  frequency  added  mass  and  moment  of  inertia.  This  section  defines  these  terms  and 
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explains  why  they  are  used. 

LAMP  calculates  the  total  instantaneous  force,  Fo,  and  moment,  Mo,  for  the  ship  rigid 
body  dynamic  equations.  Force  and  moment  contributions  from  hydrostatics,  the  incident 
wave,  hydrodynamics,  external  forces,  and  body  forces  are  used  to  calculate  Fo  and  Mo.  Due 
to  the  total  nature  of  the  force  and  moment  calculation  the  infinite  frequency  added  mass  and 
moment  of  inertia  for  the  ship  are  not  required  in  the  dynamic  equations.  However,  they  are 
used  in  the  LAMP  rigid  body  dynamic  solution  for  numerical  stability. 

The  added  mass  and  added  moment  of  inertia  terms  are  referred  to  the  global  coordinate 
system  when  calculated.  This  causes  them  to  be  time-varying  as  the  rigid  body  orientation 
changes  with  respect  to  the  global  system.  The  added  mass  and  added  moment  of  inertia  terms 
are  defined  at  an  each  time  instant  as 


aa-pfJk^dS  i,j  =  1.2 . 6  (2.26) 

B 

where  =  =  1, 2, 3  and  =  (K  x  n>)^3,  *  =  4, 5, 6.  Vector  rt  is  the  position  vector 

on  the  body,  B,  in  the  local  coordinate  system  and  ~n  the  outward  normal  of  the  fluid  on  B. 
The  global  added  mass  matrix,  Ao,  is  a  6x6  matrix  constructed  with  term  i,j  of  Ao  equal  to 
aij.  ■  The  terms  in  Ao  are  defined  as  3x3  submatrices 


Ao  = 


Ano  A120 
A210  A220 


(2.27) 


2.3.2  Translation: 

By  replacing  the  rigid  body  mass  matrix,  Mo,  with  Mo  +  A the  equation  of  motion  for 
translation  can  be  written  as, 


_ ^  J  _  _ ^ 

Fo  =  —((Mo  +  A m(t))(V0  -  Uship )) 
at 


(2.28) 
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For  stability  in  the  Runge-Kutta  numerical  integration  scheme,  added  mass  terms  are  added  to 
both  sides  of  equation  2.28  to  obtain 


Fo  +  Aiio^(Vo)  +  Ai20^(  w')  =  Jj(( Mo  +  A  m(t))(V*0  -  U  ship))  +  An  o^(v<>)  +  A12o^(v) 
1,1  at  *  (2.29) 

Rewrite  equation  2.29  by  substituting  "F  for  the  left  hand  side  where  F  =  Fo  +  Ano^(V0)  + 

Ai»&(«?) 

^F  =  ^((Mo  +  Am(t))(V0  -  17 ship))  4-  Ano— (VQ)  +  Ai2o^(^)  (2.30) 

then  group  terms  and  carry  out  differentiation  to  get  the  final  form  of  the  translation  equation 
of  motion 

F  =  (Mo  +  Am(<)  +  AMo)ft(V0)  +  A^IT|(^/)  +  (V^  -  Uship)jt(Am(t))  (2.31) 

2.3.3  Rotation: 

The  rotation  equation  of  motion  is  expressed  in  the  local  frame  Or.  The  local  frame  is  used  so 
that  derivatives  of  I  and  L  do  not  have  to  be  calculated. 

Definitions  and  relations: 

Some  terms  need  definition  prior  to  developing  the  rotation  equation  of  motion.  I  is  the  rigid 
body  mass  moment  of  inertia  tensor  in  the  Og  or  Og/  frame.  H  is  the  angular  momentum 
about  eg  in  the  Og  or  Og/  frame.  H  /  is  the  angular  momentum  about  eg  in  the  Or  frame. 


The  following  equations  relate  the  rotational  terms: 

H  =  IuT 

(2.32) 

H/  =  I/a?/ 

(2.33) 

Mo/  =  LMo 

(2.34) 

I  =  XrI/X 

(2.35) 
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Formulation: 

The  equation  of  motion  for  rotation  about  eg  in  the  O/  coordinate  system  can  be  written  as, 

Mo/  =  4(H/)  (2.36) 

at 

In  order  to  expand  the  term  ^(H/),the  following  equation,  which  is  a  standard  result  from  any 
dynamics  textbook  for  vector  derivatives  in  rotating  systems,  is  used, 

|(A)  =  |(AWxA  (2.37) 

In  equation  2.37,  A)  is  the  absolute  rate  of  change  of  A  written  terms  of  the  unit  vectors 
in  the  rotating  system,  -^(A)r  is  the  rate  of  change  of  A  as  viewed  from  the  rotating  system, 
and  Tj  is  the  absolute  angular  velocity  of  the  rotating  system.  Placing  the  angular  momentum 
vector  H/  into  equation  2.37  gives, 

4(H/)  =  4(H')r  +  <7/  X  H/  (2.38) 

Ctt  Cut 

The  angular  momentum  as  viewed  from  the  rotating  (local)  system  is, 

(H/)r  =  I/cj/  (2.39) 

inserting  equations  2.33,  2.36,  and  2.39  into  2.38  gives, 

Mo/  =  4(i'^')  +  w7  x  I/w7)  (2.40) 

at 

Due  to  time  varying  mass,  the  rigid  body  mass  moment  of  inertia  tensor,!/,  is  replaced  by 
I of  +  AI/(t).  Applying  equation  2.34  to  2.40  and  making  the  substitution  for  the  time  varying 
inertia  results  in 

LMo  =  1((I0/  +  AI f(t))u}f)  +  ~u/x  (Lot  +  AV(t))Z}/)  (2.41) 
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Then,  by  performing  differentiation  and  moving  L  to  the  right  hand  side  equation  2.41  becomes 

Mo  =  LT  (T of  +  Al/(f))^cj/  +  -^(Al /(t))ur  +  Zj/x  [(To/  +  Al/(f))  w '/]  (2.42) 

Similar  to  equation  2.29,  for  stability  of  the  numerical  scheme,  added  mass  terms  are  added  to 
both  sides  of  the  equation  2.42.  Also,  the  substitution  of  M  =  Mo  +  A2io(§(V0)  +  A220^(  £;) 
is  made  so  that 


M  =  LT  (Io/  +  AI/(t))^~u>7  +  —  (AI/(f))  to  /  4-  w7  x  [(To/  +  AI/(t))  wV]  +A2i0— (v!)+A220— (A") 


(2.43) 


Finally,  grouping  terms  gives  the  final  form  of  the  rotation  equation  of  motion 


M  =  A2io^(Vo)+  [a220LT  +  LT(lof  +  Al/(t))J  —  w/+LT— (AI/(t))w‘/+LT(  w/x  [(Io/  +  AI/(f))  w  /]) 

(2.44) 


2.3.4  Coupled  Rotation  and  Translation  Equations: 

The  translational  and  rotational  equations  of  motion,  equations  2.31  and  2.44,  are  a  coupled 
system 

1  d 
E\  Tt 

where 

1  Ano  “h  Mo  +  A7n(i)  A120L 

E  ~ 

A2 io  A22 qL  +  L  (Io/  +  AI/(f)) 

and 

I  if  —  (Vo  -  Uship)—(Am(t)) 

■  Q  J  M  -  LT( uf  x  (To/  +  AT r(t))c?r))  -  LT f 

The  details  of  [E],  and  [g]  were  not  shown  in  equation  2.24.  They  contain  terms  similar  to 
equation  2.45  except  that  there  are  some  new  terms  in  equation  2.45  due  to  the  time-dependent 
mass  and  mass  moment  of  inertia.  The  new  terms  are  (V„— t/ ship)-^(Am(t))  in  the  translation 
equation  and  LT  «7^(AT/(t))  in  the  rotation  equation.  Also,  mass  and  mass  moment  of  inertia 
vary  with  time  in  equation  2.45  but  are  constant  in  2.24. 


(2.46) 
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Chapter  3 


Models  for  Flooding 


3.1  Compart  mentation 

Watertight  internal  subdivision  using  bulkheads  to  form  compartments  within  a  ship  has  been 
the  primary  means  of  limiting  the  extent  of  flooding  in  a  damaged  ship.  During  ship  design, 
various  bulkhead  arrangements  are  evaluated  against  operational  and  damaged  stability  re¬ 
quirements  to  determine  an  optimal  ship  subdivision.  In  general,  improved  damaged  stability 
performance  with  increased  subdivision  must  be  balanced  against  drawbacks  to  compartmen- 
tation  such  as  weight,  interference  with  arrangements,  and  access  to  systems. 

The  capability  to  compartmentalize  hull  geometry  had  to  be  added  to  LAMP  in  order  to 
model  flooding  for  damaged  stability  analysis.  Since  this  thesis  is  concerned  with  the  basic 
changes  necessary  to  LAMP  for  it  to  be  used  as  a  damaged  stability  prediction  tool,  it  was 
considered  adequate  that  the  compartmentation  model  only  use  transverse  bulkheads.  If 
required  for  a  specific  ship  configuration  or  damaged  scenario,  more  detailed  compartmentation 
model  with  longitudinal  bulkheads  and  damage  control  decks  could  be  added. 

A  LAMP  program  module  was  written  to  accept  arbitrary  transverse  bulkhead  locations 
as  input.  The  bulkhead  locations  are  specified  by  their  distance  from  the  ship  bow.  The 
distance  is  normalized  through  dividing  by  the  overall  hull  length.  Any  number  of  bulkheads 
can  be  created.  The  bulkhead  locations  are  then  automatically  spliced  into  the  hull  geometry 
description  to  form  compartments  that  are  bounded  by  forward  and  aft  bulkheads,  the  hull,  and 
the  weatherdeck.  Compartments  at  the  bow  or  stem  of  the  ship  use  the  hull  geometry  to  serve 
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the  place  of  a  forward  or  aft  bulkhead  as  applicable.  Figure  3-1  illustrates  compartmentation 
of  the  bow  of  a  DDG51  Class  hull.  In  this  figure  the  intersection  of  each  line  with  other  lines 


Figure  3-1:  Compartmentation  Model  of  a  DDG51  Class  Bow 

is  described  by  an  three  dimensional  coordinate  point. 

Consistency  between  the  flooded  volume  calculation  and  the  LAMP  hydrostatic  calculation 
is  important  for  an  accurate  representation  of  the  ship’s  mean  body  position.  Because  the 
compartmentation  module  uses  the  detailed  LAMP  hull  geometry  description  to  define  the 
majority  of  a  compartment,  the  model  provides  excellent  consistency  between  the  compartment 
flooded  volume  calculation  and  the  hydrostatic  calculation.  Results  were  obtained  for  a  ship’s 
final  hydrostatic  position  after  a  flooding  event.  A  comparison  of  the  flooded  volume  against 
the  resulting  change  in  ship  displaced  volume  showed  the  two  agreed  to  within  99%. 

Finally,  the  compartmentation  module  is  only  used  for  calculations  on  the  internal  flooded 
water  and  does  not  affect  the  hull  geometry  description.  This  is  important  because  the  hull 
geometry  description  defines  the  body  on  which  the  hydrodynamic  calculations  are  performed. 
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3.2  Calculations  For  a  Compartment’s  Flooded  Volume 


For  sloshing  of  a  flooded  volume  of  water  it  was  conservatively  assumed  the  water  within  a 
rolling  compartment  maintains  a  horizontal  surface.  In  actual  ship  compartments  there  are 
generally  some  solid  objects  that  will  project  through  the  surface  of  the  water  to  reduce  free 
surface  motion.  Ship  designers  account  for  this  effect  through  the  surface-permeability  factor. 

The  horizontal  surface  assumption  produces  a  first  order  approximation  for  sloshing.  For 
simplicity,  pitch  motion  was  not  considered  in  the  flooded  water  sloshing  calculations.  Damaged 
stability  criteria  mainly  look  at  transverse  stability  while  in  a  damaged  condition  and  sloshing 
due  to  pitch  would  have  little-to-no  effect  on  this  stability.  Ignoring  pitch  motion  also  made  it 
easier  to  check  the  accuracy  of  the  flooded  volume  calculations. 

Numerical  solution  techniques  for  calculating  the  instantaneous  dynamic  free  surface  in  a 
flooded  compartment  are  an  alternative  to  assuming  a  horizontal  surface  but  they  provide  a 
relatively  small  increase  in  accuracy  compared  to  the  significant  increase  in  complexity  and 
calculation  time.  An  actual  ship  compartment  is  outfitted  and  filled  with  equipment  so  even  if 
the  instantaneous  free  surface  were  calculated  it  probably  would  not  reflect  the  actual  conditions 
in  the  flooded  compartment.  Efforts  to  accurately  calculate  a  free  surface  are  better  suited  for 
the  shallow  water  flow  problem  that  arises  when  shipping  water. 

Figure  3-2  illustrates  the  instantaneous  position  of  a  flooded  compartment  volume  with  the 
ship  undergoing  a  roll  to  starboard  with  roll  angle,  $r.  The  roll  is  made  through  the  ship  s 
center  of  gravity,  eg.  The  unprimed  coordinate  system  in  the  figure,  the  Y Z  system,  is  the 
LAMP  initial  static  system  which  is  used  as  the  reference  to  describe  the  hull  and  compartment 
geometry  through  position  vector  Rxyz.  The  primed  coordinate  system  has  been  introduced 
for  the  flooded  volume  calculation  and  rotates  to  maintain  the  Yt  axis  parallel  with  the  water 
surface. 

Flooded  compartment  volume  is  calculated  in  the  Y/Z/  system  by  summing  the  rectangular 
parallelepipeds  of  length  delyf  and  height  dzt.  These  are  indicated  in  figure  3-2  and  have  a  depth 
dx  along  the  hull’s  longitudinal  axis.  The  parallelepiped  geometry  was  selected  to  simplify 
moment,  of  inertia  calculation.  The  flooded  volume  calculation  is  initiated  by  converting  the 
compartment  geometry,  R,  to  the  primed  coordinate  system,  RJ,  where  C  is  the  rotation  matrix 
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Figure  3-2:  Compartment  Flooded  Volume  Model 


between  YZ  and  YfZf , 


Ri  =  C(Rxyz  —  Rgrav) 


(3.1) 


Since  only  roll  motion  is  considered,  the  x  and  xt  values  are  the  same.  If  it  were  desired  to 
include  pitch  motion  in  the  sloshing  model,  the  pitch  angle  could  be  included  in  C  and  the 
calculation  proceed  in  the  same  manner  as  described  here.  The  summation  for  volume  is  made 
over  x  using  j  =  jmax  intervals  and  over  zf  with  k  =  fcmax  intervals  where  x  spans  the  flooded 
compartment  length  and  zt  ranges  from  the  minimum  value  to  the  waterline  wl 


Vcmpt  =  Y  delV,dz,dx 

Xj  Zlk 


(3.2) 


At  a  fixed  roll  angle,  the  Waterline  value  wl  serves  to  determine  compartment  flooded  volume. 
If  a  specific  instantaneous  flooded  volume  is  desired  for  a  particular  flooding  scenario,  wl  can 
be  adjusted  through  an  iteration  process  until  the  desired  volume  is  reached. 

The  center  of  mass  of  a  flooded  compartment,  CGcmpt  =  (x/j  ,x/2,xfs),  can  be  calculated 


Equation  3.1  is  then  used  to  refer  CGCmpt  to  the  YZ  coordinate  system.  With  the  flooded 
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volume  and.  center  of  mass  known,  forces,  moments,  and  the  time-dependent  mass  terms  can 
be  calculated  due  to  the  flooded  water  and  added  to  the  dynamic  equations  of  motion  2.45.  If 
the  flood  rate  is  specified  it  is  used  directly  as  the  |(Am(t))  term  in  the  equations  of  motion. 
Otherwise  4  (A m(t))  needs  to  be  calculated  based  on  the  flooded  volume  time  history. 

The  calculation  for  flooded  water  mass  moment  of  inertia  is  made  by  summing  the  mass 
moment  of  inertia  of  each  parallelepiped  about  the  flooded  volume  center  of  mass.  Appendix 
A  outlines  the  method  for  doing  this.  After  calculating  the  flooded  volume  mass  moment  of 
inertia  tensor  about  the  center  of  mass,  S M I cmpt,  the  parallel  axis  theorem  and  a  rotational 
transformation  must  be  applied  to  refer  SMIcmpt  to  the  ship’s  eg  in  the  local  coordinate  system. 
The  instantaneous  flooded  volume  mass  moment  of  inertia  tensor  is  used  for  the  AI f(t)  term 
in  the  dynamic  equations  of  motion.  The  ^(AI/(t))  term  must  be  calculated  based  on  the 
flooded  volume  mass  moment  of  inertia  time  history. 

To  better  model  an  actual  ship,  compartment  permeability  could  be  used  in  the  calculations 
for  a  compartment’s  flooded  volume.  It  would  be  a  trivial  matter  to  add  permeability  to  the 
calculations  above. 

3.3  Flooding  Simulation 

With  the  compartment  flooded  volume  model  established  there  are  several  approaches  to  run¬ 
ning  a  time  domain  flooding  simulation.  First,  the  simulation  can  start  with  the  initial  condition 
that  the  flooding  event  is  complete  and  the  ship  is  in  its  final  flooded  static  condition.  A  com¬ 
puter  module  was  written  to  run  this  type  of  simulation.  The  module  solves  for  a  ship’s  final 
static  position  after  a  flooding  event  in  any  specified  compartments  is  complete.  The  module 
uses  an  iterative  procedure  to  calculate  the  final  ship  position  and  then  revises  the  ship  mass, 
moment  of  inertia  tensor,  and  eg  location  to  account  for  the  added  water.  Alternatively,  the 
module  can  maintain  the  ship  intact  and  provide  the  total  water  that  would  be  added  if  specific 
compartments  were  flooded.  This  information  can  be  used  as  an  upper  bound  on  flooded 
compartment  volume  if  it  is  desired  to  flood  the  ship  as  time  advances  in  the  simulation. 

The  second  approach  for  a  time  domain  flooding  simulation  is  to  start  with  either  an  intact 
condition  or  some  fully  flooded  compartments  and  then  flood  additional  compartments  as  the 
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simulation  time  progresses.  This  approach  is  what  is  referred  to  as  progressive  flooding  in 
Chapter  1.  The  mass  addition  rate  from  flooding  would  need  to  be  specified  in  this  type  of 
simulation.  One  technique  to  specify  this  rate  is  to  specify  a  hull  opening  location  due  to 
damage  and  let  water  enter  the  compartment  when  the  instantaneous  free  surface  is  above  the 
hole.  The  flow  through  the  hole  of  area  A  and  at  static  head  h  would  be  governed  by  the 
short  tube  orifice  equation,  Q  =  CdAy/2gh.  Tables  may  be  obtained  for  the  coefficient  Cd  from 
textbooks  on  hydraulics  such  as  LeConte  [18].  In  [6]  Dillingham  states  Cd  may  be  taken  as 
0.60  with  very  good  accuracy. 


3.4  Wind 


LAMP  is  configured  to  include  external  forces  in  its  calculation  and  currently  has  modules 
that  calculate  external  forces  for  items  such  as  appendages,  viscous  roll  damping,  and  moving 
weights.  It  is  a  simple  matter  to  include  a  heeling  moment  caused  by  beam  winds  using  an 
equation  such  as  the  following  from  reference  [17] 


Mw 


A(Vfy)2AZ(  cos($r))2 


(3.4) 


In  equation  3.4  Vw  is  the  wind  velocity,  A  is  the  ship  sail  area,  l  is  level  arm  from  centroid  of  sail 
area  to  half  draft,  4>r  is  the  roll  angle,  and  A  is  the  ship  displacement.  A  is  a  constant  whose 
value  can  vary  depending  on  units  used  in  the  equation  and  on  assumptions  on  values  for  wind 
drag  coefficient.  Reference  [17]  contains  a  discussion  on  calculating  wind  heeling  moments. 
Equation  3.4  could  be  expanded  to  include  the  wind  heading  angle  so  that  three  dimensional 
wind  forces  and  moments  on  the  ship  are  included  in  the  LAMP  calculation. 


3.5  Causes  of  Loss  of  Accuracy  in  LAMP  Flooding  Simulations 

When  linear  hydrodynamics  are  used  in  LAMP  to  solve  the  potential  flow  problem  (LAMP-1  and 
LAMP-2)  the  initial  mean  body  boundary  position  is  used  for  the  duration  of  the  calculation. 
However,  due  to  sinkage  and  trim  from  flooding  the  actual  mean  body  boundary  position  will 
change.  As  the  flooding  ship’s  mean  body  position  diverges  from  the  initial  position  a  loss  of 
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accuracy  in  the  calculated  ship  motions  will  result. 

One  strategy  to  limit  the  loss  of  accuracy  in  a  flooding  simulation  using  linear  hydrodynamics 
is  to  select  an  intermediate  body  position  between  the  intact  and  final  flooded  conditions. 
Another  approach  would  be  to  run  the  linear  hydrodynamics  LAMP  simulation  at  the  intact 
body  position  and  then  at  the  final  flooded  body  position  and  use  the  worse  case  ship  motions 
as  the  motion  estimate.  Of  course,  non-linear  hydrodynamics  (LAMP-4)  could  be  used  for 
the  flooding  simulation  but  the  penalty  is  that  much  more  time  is  required  to  perform  the 
calculation. 

As  water  floods  into  the  ship  it  causes  a  time-dependent  shift  in  the  ship’s  center  of  gravity, 
eg.  This  shift  is  not  accounted  for  in  the  LAMP  dynamic  equations  of  motion  solution.  For  a 
small  amount  of  flooding  in  a  massive  ship  the  shift  in  eg  will  be  trivial.  Under  these  conditions 
the  calculated  ship  motions  would  be  reasonably  accurate.  However,  as  the  magnitude  of  the 
eg  shift  increases,  the  calculated  angular  ship  motions  will  be  wrong  because  the  rigid  body 
moment  determined  at  each  time  step  grows  in  error. 

Finally,  when  simulating  significant  hull  damage  such  as  a  compartment  size  hole,  consid¬ 
eration  should  be  given  as  to  whether  complete  body  panelization  that  assumes  an  intact  hull 
will  provide  an  accurate  enough  hydrodynamic  solution.  A  more  accurate  method  may  be  to 
re-panelize  the  hull  around  the  physical  damage. 
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Chapter  4 


Models  for  the  Green  Water 
Problem 


4.1  Background  and  Scope  of  Green  Water  Model 

As  a  flooding  ship  loses  freeboard,  the  likelihood  of  shipping  water  increases.  The  water  on 
deck,  in  turn,  can  further  degrade  the  ships  stability  condition.  This  is  of  special  concern  for 
smaller  vessels.  A  review  paper  on  the  subject  of  water  on  deck  and  stability  of  has  been 
published  in  reference  [3].  Also,  in  reference  [6],  Dillingham  provides  calculations  to  show  that 
instabilities  caused  by  excessive  deck  water  may  cause  capsizing  of  small  fishing  vessels.  The 
green  water  problem  can  be  decomposed  into  three  subproblems:  shipping  water,  water  escaping 
off  deck,  and  the  motion  of  water  on  deck.  For  adequate  modeling  of  the  green  water  problem 
in  a  time  domain  ship  motion  program  each  subproblem  must  be  addressed. 

Most  work  reported  on  the  green  water  problem  has  focused  on  the  water  motion  on  deck 
subproblem.  This  emphasis  over  the  other  two  subproblems  is  primarily  due  to  water-motion- 
on-deck  sharing  the  same  basic  theory  as  that  for  flow  of  compressible  gases.  Specifically, 
the  equations  governing  water  motion  on  deck,  equation  4.13,  are  derived  from  the  theory  for 
waves  in  shallow  water,  covered  in  detail  in  reference  [29],  and  are  of  the  same  form  as  the 
compressible  gas  dynamics  equations.  The  necessity  in  industry  for  dealing  with  the  flow  of 
gases  has  resulted  in  numerical  solution  methods  that  can  be  directly  applied  to  the  water- 
motion-on-deck  problem.  The  objective  in  solving  the  shallow  water  equations  is  to  solve  for 
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the  water  depth  and  horizontal  particle  velocity.  Solutions  to  the  governing  equations  can 
involve  discontinuities  such  as  shocks,  bores,  and  hydraulic  jumps.  Any  numerical  methods 
used  to  solve  the  equations  must  be  capable  of  treating  discontinuities.  Most  of  the  numerical 
methods  fall  into  two  general  categories  characterized  by  the  scheme  used  to  solve  the  governing 
equations:  the  flux  difference  splitting  method  and  the  random  choice  method,  also  known  as 
Glimms  Method.  Both  schemes  handle  discontinuities  in  the  solution  without  any  special 
treatment  and  were  evaluated  in  this  thesis.  The  Random  Choice  method  involves  looking 
at  solution  curves,  called  characteristics,  to  the  governing  differential  equations.  The  flux 
difference  splitting  method  combines  a  finite  difference  method  with  characteristics.  A  different 
approach  to  solving  the  water  on  deck  problem  is  in  reference  [13]  where  the  authors  solve  for 
wave  motion  in  a  rolling  tank  using  a  finite  difference  scheme  coupled  with  analytical  techniques. 
This  approach  illustrates  the  special  care  that  must  be  taken  with  discontinuities  with  which 
the  numerical  scheme  is  incapable  of  handling. 

This  chapter  develops  a  green  water  model  for  the  LAMP  program  in  head  seas.  The  two 
main  methods  for  solving  the  water-on-deck  subproblem  are  evaluated  for  use  in  the  model  using 
a  two-dimensional  free  surface.  A  method  for  incorporating  water  shipping  into  LAMP  is  also 
devised.  The  water  escaping  off  deck  subproblem  is  not  formulated  because  a  proper  calculation 
would  require  a  three-dimensional  free  surface  solution  to  the  shallow- water  problem.  Three 
dimensions  would  introduce  transverse  water  velocities  so  that  as  the  green  water  travels  aft  on 
the  weatherdeck  it  also  moves  towards  the  port  and  starboard  deck  edges  and  then  over  board. 
The  green  water  model  used  for  this  thesis  only  calculates  longitudinal  water  velocities  so  green 
water  mass  is  removed  from  the  weatherdeck  by  letting  it  fall  off  the  after  end  of  the  portion 
of  the  weather  deck  included  in  the  computation.  This  is  accomplished  by  setting  boundary 
conditions  for  the  water-on-deck  calculation  to  zero  at  the  aft  end  of  the  weatherdeck 

4.2  Flux  Difference  Splitting  Method  for  Water  Motion  on  Deck 

The  flux  difference  splitting  method  was  developed  based  on  the  flux  vector  splitting  method 
originally  introduced  by  Steger  and  Warming  in  [28].  Steger  and  Warming  developed  a  basic 
theory  of  the  flux  vector  splitting  method  to  compute  the  shock  wave  for  the  gas  dynamic  equa- 
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tions.  Because  of  differences  between  the  non-homogeneous  governing  equations  for  shallow 
water  flow  on  deck  and  gas  dynamics,  it  is  better  when  solving  the  shallow  water  equations 
to  split  flux  differences  instead  of  the  flux  vector.  In  reference  [l]  Alcrudo  presents  the  flux 
difference  splitting  method  to  solve  problems  for  open  channel  hydraulics.  In  references  [12] 
and  [11]  the  flux  difference  splitting  method  is  applied  to  solve  the  shallow  water  flow  on  deck 
problem. 

The  flux  difference  splitting  method  is  a  an  upwind  (one-sided)  finite  difference  scheme  that 
solves  the  shallow  water  wave  propagation  problem  for  a  two  dimensional  free  surface.  For 
supercritical  flow,  wave  information  can  only  travel  downstream.  In  the  case  of  subcritical 
flow,  wave  information  will  travel  in  both  directions.  In  order  to  construct  an  upwind  scheme 
valid  for  all  regimes  and  directions  of  flow,  a  decomposition  of  the  flux  related  to  positive 
and  negative  propagation  speeds  is  needed.  The  flux  decomposition  is  devised  so  that  wave 
information  can  not  travel  upstream  in  supercritical  flow.  For  flux  difference  splitting  methods 
the  flux  difference  operator,  A~F  and  ARHS  in  the  equations  below,  is  split  based  on  the 
characteristic  directions. 

Reference  [12]  can  be  used  to  formulate  the  shallow  water  flow  governing  equations  for  a 
two-dimensional  free  surface  with  the  geometry  illustrated  in  Figure  4-1.  Many  of  the  variables 
in  the  governing  equations  are  not  defined  in  the  figure.  Of  these,  u\  is  the  ship  surge  velocity, 
u3  is  the  ship  heave  velocity,  ©  is  the  pitch  angle,  u5  is  the  angular  velocity,  and  g  is  acceleration 
due  to  gravity.  Note  that  u  is  the  water  particle  velocity  in  the  x-direction.  The  governing 
equations  in  vector  form  is  _  __ 


where  W  = 


and  C  =  I 

y  94PC  ) 
They  are  defined  as 


ugC 

u2gC  +  ^(gC)2 


h(gC)2 


■\,[D}  = 


o  qigC  +  Q2U  +  93 


.  The  variables  q\  through  q\  are  functions  of  deck  motion  and  geometry. 
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Figure  4-1:  Coordinate  System 


for  Two-Dimensional  Free  Surface 
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Equation  4.1  is  derived  by  applying  the  shallow  water  assumptions  to  the  continuity  equation 
and  Eulers  equations  of  motion. 

The  derivatives  of  the  flux  vector  ~F  and  #  can  be  expressed  in  terms  of  W 


8~F  rj.dW 
dx  ~  ^  dx 


and 


dH  '7,dW 
dx  ~  *  ^  dx 


(4.2) 


where  [Ji]  and  [J2]  are  the  Jacobian  matrices.  The  eigenvalues  of  [Ji]  are 


Xi  =  u  +  y/gC  and  \2  =  u-  1/9C 


(4.3) 


with  the  eigenvectors 


and 


(4.4) 
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The  difference  of  W  and  flux  vector  ~F  are  approximated  as 

2  2 

Atf^  =  Yaket  and  A  F  =  A  kQ-k^k  (4-5) 

k= 1  k= 1 

The  right  hand  side  of  equation  4.1  equal  to  zero  corresponds  to  water  sloshing  for  a  given 
initial  water  surface  profile  and  the  ship  stationary.  With  the  right  hand  side  zero,  the  finite 
difference  scheme  with  the  split  flux  difference  can  be  expressed  in  the  following  form 


where 

F*+±  =  2 J+1  +  F i)  ~  2  ^  ak>iH  K-i+i  (4-7^ 

F;_1  =  i  +  F3~ i)  “  2  |\i-|  (4-8) 

The  scheme  in  equation  4.6  is  of  the  first  order. 

When  the  right  hand  side  of  equation  4.1  is  not  equal  to  zero,  it  is  treated  as  another  flux 
difference  term  and  projected  into  the  eigenvector  space  as  follows 

oTt  i  2 

is3=[£l|^r  +  c  =  ~A (49) 

The  right  hand  side  flux  difference,  where  the  flux  difference  is  related  to  equation  4.9  by 
A  RHS  =  RFISAx,  is  then  split 

“d  (41°) 


where 


xij-i  =  and  \j+i~  2^+5  I Afc.j+?|)  ^■11^ 


The  spht  flux  difference  is  then  included  in  the  finite  difference  scheme  of  equation  4.6 

^n+,) = wf + #(^-i  +  aSS2;+p  (4-12) 
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The  minus  sign  in  equation  4.12  prior  to  the  R.H  S  term  is  due  to  the  minus  sign  in  equation 
4.9.  In  [28]  Steger  shows  that  flux  splitting  schemes  are  stable  if  and  only  if  A*  <  1. 

Shallow  water  first-order  schemes  suffer  from  numerical  dissipation  (the  shock  front  will  be 
smeared).  Reference  [12]  shows  how  a  flux  limiter,  which  is  a  correction  term  for  the  numerical 
flux,  can  be  applied  to  the  first  order  scheme  to  make  the  scheme  higher  order. 

The  flux  difference  splitting  method  can  be  used  to  solve  the  shallow  water  equations  for  a 
three-dimensional  free  surface.  Reference  [12]  formulates  a  technique  using  the  flux  difference 
splitting  method,  together  with  the  Fractional  Step  Method  [32],  so  that  solutions  to  the  shallow 
water  equations  can  be  obtained  by  solving  two  sets  of  two-dimensional  free  surface  problems. 
The  three-dimensional  solution  works  as  follows.  Two  sets  of  vector  equations  similar  to 
equation  4.1  are  devised.  One  set  is  for  flow  in  the  ship  longitudinal  direction  (x  -  axis ) 
where  each  specific  equation  holds  for  a  certain  deck  strip  of  thickness  Ay.  The  other  set  of 
vector  equations  is  for  flow  in  the  transverse  direction  (  y  axis )  where  each  equation  is  for  a 
certain  deck  strip  of  thickness  Ax.  Then  the  Fractional  Step  Method  advances  the  solution  in 
time.  At  each  time  step  the  sets  of  equations  along  the  x  -  axis  are  solved  for  an  intermediate 
solution  assuming  no  y  dependency.  Then,  in  the  same  time  step,  the  sets  of  equations  along 
the  y  -  axis  are  solved  from  the  intermediate  solution  assuming  no  x  dependency.  Instead 
of  solving  the  three-dimensional  free  surface  governing  equation  on  (m  x  n)  nodes,  a  total  of 
(m+n)  two-dimensional  free  surface  equations  are  solved  along  the  x  and  y  directions  separately 
using  the  Fractional  Step  Method. 

The  difference  schemes  of  equations  4.6  and  4.12  for  a  two-dimensional  free  surface  were 
programmed  in  a  computer  so  that  the  method  could  be  compared  with  Glimms  method. 

4.3  Glimms  Method  (Random  Choice  Method)  for  Water  Mo¬ 
tion  on  Deck 

For  a  two-dimensional  free  surface,  Glimms  Method  is  performed  by  dividing  the  physical 
domain  into  intervals,  i  —  1,2, 3, ...,  i  max.  In  each  interval  at  time  nA t  the  solution  is  approxi¬ 
mated  by  piecewise  constant  depth,  Ci;  &nd  particle  velocity,  uz .  At  the  boundary  between  each 
interval,  the  depth  and  velocity  are  therefore  discontinuous  which  gives  a  sequence  of  Riemann 
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problems  which  may  be  solved  to  advance  the  solution  to  time  (n  +  l)At.  The  solution  to 
the  Riemann  problem  is  outlined  below  and  the  full  solution  is  developed  in  Appendix  B.  The 
Riemann  problem  is  also  known  as  the  dam  breaking  problem.  Ftojn  the  time  (n  -1-  l)At  Rie¬ 
mann  solutions,  it  is  desired  to  construct  another  piecewise  solution.  This  is  done  by  randomly 
sampling  the  solution  in  each  interval  and  then  using  the  sampled  value  as  an  approximation  for 
that  interval’s  constant  piecewise  solution.  In  [10]  Glimm  showed  the  random  sampling  scheme 
converges  to  a  weak  solution.  The  disturbance  resulting  from  the  solution  for  each  individual 
Riemann  Problem  must  not  be  able  overlap  with  the  disturbance  from  the  adjacent  Riemann 
Problem.  No  disturbance  must  be  allowed  to  propagate  further  than  Amt*rval  in  a  time  At. 
Thus  the  resulting  Courant  condition  that  must  be  satisfied  is  A*2A T3^  •>  (M  "b  V^C)-  The 
solution  obtained  is  unconditionally  stable. 

Glimms  Method  is  used  to  solve  the  water  motion  on  deck  problem  in  references  [6],  [22], 
and  [23].  References  [27]  and  [4]  illustrate  use  of  the  Glimms  Method  to  solve  gas  dynamic 
problems. 

Appendix  C  contains  a  MATLAB  program  that  solves  a  Riemann  Problem  using  Glimms 
Method. 


4.3.1  Solution  of  the  Riemann  Problem 


Glimms  method  requires  solution  of  the  Riemann  Problem,  also  known  as  the  dam  breaking 
problem,  at  each  spatial  interval  in  the  computational  domain  at  each  time  interval.  A  summary 
of  the  Riemann  Problem  solution  procedure,  taken  from  Dillingham  [6]  follows.  Appendix  B 
provides  details  of  the  calculations.  In  keeping  with  Dillingham’s  notation,  the  y  axis  is  used 
as  the  coordinate  axis  for  the  Riemann  Problem,  v  is  the  horizontal  velocity,  and  A  is  used  for 


the  water  depth. 

The  Riemann  Problem,  illustrated  in  figure  4-2,  consists  of  solving  a  system  of  nonlinear 
hyperbolic  equations,  called  the  shallow  water  wave  equations, 


dv  dv  dX 

~xt  4-  v~x — b  9xt-  —  0 
dt  dy  dy 


8X  dX  dv 

a"d  m+vTy+xa iT° 


(4.13) 


for  a  two-dimensional  free  surface  where  g  is  gravity  acceleration,  subject  to  initial  conditions 
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on  the  right-hand  and  left-hand  sides 


v  =  V0 
A  =  A0 


(4.14) 


v  =  V<i 
A  =  A2 


)y<  0 


and 


)y>  0 


Figure  4-2:  Initial  Conditions  for  the  Riemann  Problem 


It  is  assumed  that  the  water  is  initially  higher  on  the  left  side  of  the  dam.  The  dam  is  then 
removed  at  time  t  =  0.  Depending  on  initial  conditions,  the  Riemann  problem  solution  falls  into 
different  categories,  identified  as  Cases  I  through  V.  The  variable  c  as  used  below  represents  the 
wave  propagation  speed;  it  is  the  local  speed  of  propagation  of  small  disturbances  relative 
to  the  moving  stream.  In  shallow  water  theory,  c  is  related  to  the  water  height  by  c  =  y/gX. 


Solution  to  the  Riemann  Problem:  Case  I 
If  V2  +  2c2  >  Vo  +  2cq  and 


g  (Aq  +  Ai)(Ai  —  Ap)2  5 
2  AiAq 


(4.15) 


then  the  solution  consists  of  a  single  shock  and  a  single  rarefaction  as  indicated  in  figure  4-3Let 
co  =  VgXo,  ci  =  v^.and  c2  =  VgXi-  Also,  let  R  =  where  £  is  the  shock  speed,  then 


I 


1 

1 

V2 

V3 

1 

1 

A2 

!V1 

Vo 

A3 

1X1 

1 

_ 1 _ 

Ao 

Figure  4-3:  Solution  for  Riemann  Problem  Case  I 


solve  the  following  equation  for  R 


i 

R-Jjil  1  +  \/8^+T]  +  2  [i  (v/8i?2  +  l  -  l)] 2  = 


2  V2  -  Vo  +  2 c2 


Co 


then  we  have 


£  = 
Vi  = 

ci  = 


cqR  +  Vo 

Co  (i  +Vs^2  +  i)  +v0 

-i  l 

Co  [VSR2  +  1  -  l)  2  and  *i  =  J 


(4.16) 


(4.17) 

(4.18) 

(4.19) 


In  zone  3  we  have 

v=2^1__Vo  +  ¥l^¥°+c2'j+V0  and  A  =  ^(f2  +  2c2-|)2  (4.20) 

where  zone  3  is  bounded  by  (V2  —  c2)t  <  y  <  {Vi  —  c\) t. 

Solution  to  the  Riemann  Problem:  Case  II 

If  equation  4.15  is  not  satisfied  then  the  solution  consists  of  two  shocks  as  indicated  in  figure 
4-4 
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Figure  4-4:  Solution  for  Riemann  Problem  Case  II 


Solve  the  following  equation  for  Ai 


then 


g  (Ao  4*  Ai)(Ai  - 

-Ao)21 

1 

2 

g  (Ai  +  A2XA2  —  A1)2 

2  AiAo 

2  A2Ai 

(4.21) 


Vi 


£1 


T,  ,  J"^  (^o  +  Ai)(Ai  —  Ao)2l 2 

Vo  +  [2 - XJTo  _ 

A,  Vi  —  A(|V(|  ,  ,  A2V'2-A,Vi 

A.-A0-  a"d  A,-” 


(4.22) 

(4.23) 


Solution  to  the  Riemann  Problem:  Case  III 

If  V"0  -  P2  >2|co-C2|  then  the  solution  consists  of  two  rarefaction  waves  as  in  figure  4-5The 
solution  for  zone  1  is 

and  A1=«,(^  +  ^)2.  (4.24) 

and  A  =  -j  (v2  +  2c2  -  |)2  (4.25) 


Vi 


V'O  +  v2 


+  C2  —  Cq 


In  zone  3  we  have 


_  2  fy  K 
V  3  V*  +  2 


+  c2 
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A2 


Vo 

Ao 


Figure  4-5:  Solution  for  Riemann  Problem  Case  III 
where  zone  3  is  bounded  by  (V2  ~~  °2)  t  <  V  <  (^1  ~  ci)  I*1  Zone  4 

HIM-*)  and  A=^(2co"Vo+f' 

where  zone  4  is  bounded  by  (Vj  +  c\)t  <y  <  (Vo  +  co)  t. 


and  A  =  ^- f2co  -  Vo  +  7 
9  g  \  t, 


Solution  to  the  Riemann  Problem:  Case  IV 

If  A0  =  0  then  the  solution  is  a  single  rarefaction  wave  as  in  figure  4-6. 


Figure  4-6:  Solution  for  Riemann  Problem  Case  IV 
The  solution  in  Zone  3  is 

7  2A/  V2  \  and  A  =  If2c,+V',-^ 


where  zone  3  is  bounded  by  (V2  —  C2)  t  <  y  <  (V2  +  2C2)  t. 

Solution  to  the  Riemann  Problem:  Case  V 

If  in  case  III  there  is  the  additional  condition  that  V2  +  2c2  <  Vo  —  2co  then  the  water  depth 
will  be  equal  to  zero  in  Zone  1  and  the  problem  must  be  treated  similarly  to  Case  IV . 

4.4  Selection  of  Water  Motion  on  Deck  Method 

Computer  programs  for  the  flux  difference  splitting  method  and  for  Glimms  method  were  gen¬ 
erated  to  solve  a  Riemann  problem.  The  goal  of  this  effort  was  to  compare  the  methods  to 
see  which  is  most  suitable  for  incorporation  into  the  LAMP  water  on  deck  model.  Figure  4-7 
illustrates  the  solution  to  a  Riemann  problem  using  a  first  order  flux  difference  splitting  method 
scheme  without  flux  limiters.  The  solid  line  is  the  exact  solution,  the  other  two  lines  are  the 
numerical  solutions.  The  line  with  the  shorter  dashes  is  for  $  =  @  and  the  line  with  the 
longer  dashes  is  for  -gjj  =  The  scheme  converges  to  the  exact  solution  as  A y  is  made 

smaller. 


Figure  4-7:  Riemann  Problem  Solution  Using  Flux  Difference  Splitting  Method 

Figure  4-8  illustrates  the  solution  to  the  same  Riemann  problem  using  the  Random  choice 
method.  The  solid  line  is  the  exact  solution,  the  other  two  lines  are  the  numerical  solutions. 
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The  line  with  the  shorter  dashes  is  for  ^  and  the  line  with  the  longer  dashes  is  for 

A£  =  Mgl.  This  scheme  also  converges  to  the  exact  solution  as  A y  is  decreased. 

At/  0.01 

Due  to  the  random  interval  sampling  required  in  this  procedure,  the  results  vary  each  time 
the  calculation  is  run.  For  example,  for  constant  A t,  Ay,  and  initial  conditions,  a  set  of 
calculations  will  produce  a  unique  solution  curve  for  each  calculation.  Figure  4-9  shows  three 
different  solutions  to  a  Riemann  problem,  the  dashed  lines,  using  the  Random  Choice  Method. 
For  each  solution  the  initial  conditions  and  calculation  parameters  were  the  same. 

Other  Riemann  problems  were  investigated  to  compare  the  Random  Choice  and  first  order 
flux  difference  splitting  methods.  Table  4.1  summarizes  a  comparison  of  the  two  methods. 
An  ”X”  entered  in  the  Table  means  that  method  is  superior  to  the  other  method  in  the  given 
category. _ _ _ 


Comparison  Category 

Random  Choice  Method 

Flux  Splitting  Method 

Ease  of  Programming 

X 

Accuracy  for  Given  Discretization 

X 

Repeatability  of  Solution 

X 

Ability  to  Capture  Discontinuities 

X 

Table  4.1:  Comparison  of  Shallow  Water  Solution  Methods 


The  flux  difference  splitting  method  was  selected,  primarily  due  to  ease  of  programming,  as 
the  numerical  method  to  be  used  for  the  for  LAMP  green  water  model. 

4.5  Water  Shipping  Model 

The  difficulties  in  creating  a  realistic  water  shipping  model  are  best  framed  by  Dillingham  in 
[6].  He  states,  ’In  general  the  flow  over  the  bulwark  is  very  complicated  since  it  results  from 
the  unsteady  interaction  between  shallow  water  waves  on  the  deck  and  deep  water  waves  off 
the  deck.  In  addition,  direct  observations  of  models  indicate  that  the  deck  water  may  regularly 
impinge  on  the  bulwark  with  considerable  velocity  and  be  thrown  over  the  side  in  what  amounts 
to  a  spray.  To  refer  to  this  process  as  either  turbulent  or  nonlinear  is  an  understatement.” 

For  water  shipping  to  occur,  the  height  of  the  free  surface  at  the  ship  side  must  exceed  the 
height  of  the  bulwark  (or  deck  edge  if  there  is  no  bulwark)  and  the  relative  velocity  of  the  water 
above  the  top  of  the  deck  edge  must  be  directed  inward  onto  the  deck.  The  height  of  the  free 
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0.65 


Figure  4-8:  Riemann  Problem  Solution  Using  Random  Choice  Method 


Figure  4-9:  Solution  Randomness  in  Random  Choice  Method 


44 


surface  above  the  bulwark  or  deck  edge  is  referred  to  as  the  relative  elevation,  t]r.  This  thesis 
makes  use  of  a  ship  model  with  no  bulwark  so  hereafter  only  the  deck  edge  will  be  mentioned 
in  the  water  shipping  discussion.  Figure  4-10  illustrates  some  geometries  and  variables  used 
in  the  water  shipping  model.  The  XZ  coordinate  system  represents  the  global  system  which 
maintains  a  fixed  orientation.  The  primed  system,  XfZt ,  is  the  local  coordinate  system  that 
moves  with  the  ship.  The  variable  rj  represents  the  free  surface  elevation  and  Zd  represents  the 
height  of  the  deck  edge.  Each  is  referred  to  the  global  system  as  a  value  along  the  z  axis.  The 
instantaneous  relative  elevation  would  be  rj R(t)  =  r/(t)  —  Zd(t).  The  vectors  VWx  and  V Wy 
are  components  of  the  water  particle  velocity  V  w.  The  vectors  V  dx  and  V  dy  8X6  components 
of  the  deck  edge  velocity  ~V  d-  All  vectors  are  referred  to  the  global  system. 


Figure  4-10:  Geometry  and  Variables  for  the  Water  Shipping  Model 


4.5.1  Free  Surface  Elevation  for  Water  Shipping 

A  basic  description  of  the  instantaneous  free  surface,  rj(t),  can  be  made  by  considering  only  the 
incident  wave  potential.  However,  for  a  more  accurate  description  of  the  free  surface  around  the 
ship  hull  the  complete  incident,  diffracted,  forward  speed,  and  radiated  wave  system  potentials 
should  be  used.  The  hydrodynamic  potentials  cause  an  increase  in  free  surface  elevation  which 
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is  commonly  referred  to  as  "pile-up”  in  the  wave-body  interaction.  For  a  more  complete 
estimation  of  pile-up  the  jet  motion  of  the  free  surface  from  water-entry  and  slamming  should 
be  considered.  In  reference  [8]  a  numerical  technique  for  predicting  the  occurrence  of  water 
shipping  is  presented  and  it  is  concluded  that  the  radiation  and  diffraction  terms  play  important 
roles  in  the  water  shipping  analysis  and  cannot  be  neglected. 

For  the  water  shipping  model  presented  here,  the  incident  wave  potential  is  used  for  calcu¬ 
lating  the  free  surface  height  at  the  ship  side.  An  elevation  correction  is  then  made  to  account 
for  the  forward  speed,  radiation,  and  diffraction  potentials.  The  hull  surface  hydrodynamic 
pressure  at  the  body  panel  adjacent  to  the  waterline  and  nearest  the  desired  deck  location  is 
used  to  compute  an  elevation  which  is  added  to  the  relative  elevation  computed  from  the  ship 
motion  and  incident  wave  potential. 

For  a  LAMP-2  calculation  of  the  CG-47  hull  in  storm  seas  case,  this  elevation  correction  did 
not  prove  to  be  particularly  large  most  of  the  time.  It  is  probable  that  the  correction  would 
be  larger  for  a  LAMP-4  calculation,  as  the  hydrodynamic  calculation  would  include  some  hull 
flare  effects. 

4.5.2  Relative  Velocity  for  Water  Shipping 

Several  terms  must  be  considered  to  determine  the  relative  velocity,  V  T,  of  the  water  above  the 
deck  edge.  The  water  particle  velocity,  Vw,  and  ship  velocity  at  the  deck  edge,  V  d,  must  be 
known.  Also,  an  additional  velocity  due  to  the  effect  of  hydraulic  flow,  Vh,  must  be  included. 

Hydraulic  flow  occurs  due  to  energy  conservation  from  converting  the  static  head  from  the 
relative  elevation  into  a  velocity  according  to  Bernoulli’s  equation.  Referring  to  figure  4-11, 
where  the  variable  zf  is  the  local  height  off  the  weatherdeck,  the  static  head  height  is  (  tj  Zd )  zf 
and  the  hydraulic  velocity  is  Vh(zt)  =  ^2g(rjR  -  zf).  The  hydraulic  velocity  is  assumed  to  be 
directed  onto  the  ship  parallel  with  the  xt  axis.  The  dependency  that  Vh  has  on  z/  can  be 
removed  by  averaging  Vh  over  the  static  head.  This  is  done  with  the  following  integral  to 
obtain  mass  flow  rate,  Q, 

Q  =  d  y/lgiVR  -  z,)dz!  (4-28) 

In  [6]  Dillingham  uses  a  technique  similar  to  equation  4.28  to  integrate  for  the  mass  flow  rate 
onto  a  ship’s  weatherdeck.  Carrying  out  the  integration  in  4.28  and  then  dividing  by  the  static 
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head  height  to  obtain  an  average  hydraulic  flow,  Vhave,  gives 


ave 


(4.29) 


The  water  shipping  model,  then,  uses  Vhave,  Vw,  and  V d  for  the  relative  velocity  calculation. 

The  calculation  for  Vhave  is  a  simplification  made  for  this  thesis.  In  reference  [11]  the 
hydraulic  flow  contribution  to  water  shipping  is  treated  within  a  more  general  calculation  for 
mass  of  the  wave  shipped  on  deck.  The  mass  calculation  is  performed  by  integrating  around 
the  edge,  L,  of  the  weatherdeck  as  follows, 


Mass  =  p  r  f  O,  [jPV.  -  Vh(zl)dzl 


dldt 


(4.30) 


The  constant  CQ  in  equation  4.30  was  determined  experimentally  by  Grochowalski  to  be  0.55. 


Figure  4-11:  Water  Shipping 

Because  ~VW  and  are  vector  quantities,  the  calculation  for  V  r  can  be  performed  several 
different  ways  depending  on  the  assumptions  made.  These  calculations  will  be  referred  to  as 
Methods  I,  II,  and  III.  In  any  case,  since  shallow  water  theory  assumes  all  fluid  velocity  is 
parallel  with  the  deck  surface,  the  direction  for  V  T,  shown  in  figure  4-11,  is  always  considered 
to  be  parallel  with  the  weatherdeck.  For  Method  I,  the  calculation  will  use  only  the  horizontal 
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components  of  the  wave  and  particle  velocities  so  that 

Vr  =  Vdx-VWx+Vhave  (4.31) 

For  Method  II,  the  calculation  will  consider  the  ship’s  instantaneous  pitch  angle,  0,  and  use 
the  components  of  the  vectors  V  Wx,  V  Wy,  V  dx,  and  V  dy  that  are  parallel  with  the  x/  axis. 
The  calculation  for  Vr  in  this  case  is 

Vr  =  Vdx  cos(0)  -  Vdy  sin(0)  -  (VWx  cos(©)  -  VWy  sin(©))  +  Vhave  (4.32) 

The  Method  HI  calculation  will  be  performed  by  conserving  the  momentum  of  the  water  particle. 
For  long  crested  head  seas  there  are  only  two  wave  particle  velocity  components  so  that 

Vr  =  Vdx  cos(0)  -  Vdy  sin(0)  ±  yjvfx  +  V%y  +  Vhave  (4.33) 

where  the  plus  (+)  is  used  when  (VWx  cos(0)  -  VWy  sin(©))  <  0  and  the  minus  (-)  is  used  when 
(VWx  cos(0)  -  VWy  sin(©))  >  0. 

The  relative  velocity  and  elevation  must  be  calculated  at  each  time  step.  As  indicated  in 
figure  4-11,  the  velocity  is  considered  to  be  a  constant  from  the  deck  edge  to  the  top  of  the  free 
surface.  With  the  relative  velocity  and  elevation  of  the  water  at  the  deck  edge  known,  they 
can  be  considered  as  boundary  conditions  for  a  water  on  deck  problem.  When  the  relative 
elevation  becomes  negative,  the  boundary  conditions  are  set  to  zero.  Evolving  the  water- 
on-deck  problem  through  time,  with  the  relative  velocity  and  elevation  boundary  conditions 
updated  at  each  time  step,  will  result  in  mass  flow  onto  the  deck.  Wa  ter  shipping,  therefore, 
is  treated  as  a  boundary  condition  to  the  water  on  deck  problem. 

4.6  Green  Water  Model  in  LAMP 

To  create  a  green  water  model  in  LAMP  for  head  seas,  a  weatherdeck  was  included  in  the  hull 
geometry  description.  The  deck  was  modeled  as  a  surface  projecting  aft  from  the  ship  bow 
a  distance  approximately  25%  of  the  overall  ship  length.  The  weather  deck  is  divided  into 
strips  of  width  dely.  At  the  center  of  each  strip  is  a  deck  edge  reference  point  at  which  data 
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needed  for  the  water  shipping  problem  is  calculated.  Figure  4-12  illustrates  this  setup.  The 


dely 


Figure  4-12:  Weatherdeck  Division  for  LAMP  Green  Water  Model 

model  uses  the  flux  difference  splitting  method  to  solve  the  shallow  water  problem  in  each  strip 
as  a  two-dimensional  free  surface.  Transverse  water  flow  (from  one  strip  to  another)  is  not 
calculated. 

A  LAMP  subroutine  was  provided  by  SAIC  to  calculate,  at  each  time  step,  the  deck  edge 
motion,  deck  edge  location,  free  surface  height,  and  incident  wave  particle  velocity  for  each 
reference  point  shown  in  figure  4-12.  With  this  data,  the  relative  elevation  and  velocity  can  be 
determined  for  each  deck  edge  point. 

After  solution  of  the  shallow  water  problem  for  each  strip  at  each  time  interval,  the  forces, 
moments,  and  the  time-dependent  mass  and  mass  moment  of  inertia  terms  in  the  dynamic 
equations  of  motion  can  be  calculated  that  result  from  the  water  on  deck.  Because  of  the  time- 
dependent  mass  terms  in  the  dynamics  equations  of  motion,  the  green  water  on  deck  is  treated 
as  part  of  the  rigid  body.  Thus  the  only  green  water  induced  forces  and  moments  that  need  to 
be  calculated  for  rigid  body  motion  are  those  due  to  gravity  acceleration.  For  the  local  loads 
on  the  weatherdeck  and  superstructure,  however,  the  absolute  green  water  fluid  acceleration  in 
the  local  coordinate  system  must  be  calculated.  Let  slqw  be  the  absolute  acceleration  of  the 
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green  water  in  the  local  ship  coordinate  system.  As  defined  previously  V  is  the  rigid  body 
velocity  at  eg  and  u t  the  rigid  body  rotation.  Then  from  simple  dynamics 

a^w  =  ——  +  -7-  x"r/+  oj  x  (w  x  TV)  +2w  x  V/4-a^  (4.34) 

dt  dt 

where  "rV  and  vV  is  the  position  and  velocity,  respectively,  of  the  green  water  referred  to  the 
local  coordinate  system.  The  vector  a ^  is  the  acceleration  of  a  point  in  the  green  water 
measured  relative  to  the  local  frame  by  an  observer  attached  to  the  local  frame.  This  quantity 
is  zero  for  the  green  water  problem.  The  vertical  acceleration  component  of  a gw  is  defined  as 
a zf.  At  each  time  step  a J  can  be  calculated  for  the  green  water  on  deck  as  a  function  of  the 
x  axis  in  figure  4-12.  The  total  deck  pressure,  Pdecki  is  then  be  calculated  by  including  the 
gravity  static  head,  which  must  include  pitch  angle,  and  the  depth  of  the  green  water,  C ,  so 
that 

Pdeck  =  PC  (a*'  +  9 ooe(0))  (4-35) 

Forces  from  green  water  on  large  superstructure,  where  the  water  particle  velocity  at  the  struc¬ 
ture  boundary  can  be  approximated  as  zero  over  a  wide  area,  can  be  calculated  from  equation 
4.35  by  integrating  the  pressure  through  the  water  column.  For  smaller  weatherdeck  structures 
a  better  estimate  of  the  local  green  water  loads  may  be  made  by  considering  the  green  water 
particle  velocity. 

Figure  4-13  visualizes  a  solution  to  the  green  water  problem  by  plotting  the  water  depth  for 
one  strip  over  a  series  of  time  steps  with  a  sinusoidal  incident  wave.  In  this  case  the  ship  is 
moving  to  the  right  into  the  sea,  the  bow  is  on  the  right,  and  time  is  increasing  with  the  higher 
curves.  The  horizontal  portion  of  each  curve  is  an  area  where  the  deck  is  dry. 
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Figure  4-13:  Flow  Visualization  of  Green  Water  on  Deck 
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Chapter  5 


Validation  and  Results 


5.1  Validation 

5.1.1  Validation  of  Dynamic  Equation  of  Motion  Solver 

The  changes  to  the  LAMP  dynamic  equations  of  motion  to  include  time-dependent  mass  and 
moment  of  inertia  were  validated  using  the  model  shown  in  Figure  5-1.  The  model  is  a  body 
with  two  lumped  masses,  M0,  at  the  ends  of  a  massless  rod.  The  body  has  initial  linear  velocity 


Figure  5-1:  Model  Used  to  Validate  Dynamic  Equations  of  Motion  Solver 
VQ  and  initial  angular  velocity  u>  whose  vector  representation  is  pointed  along  the  x  axis.  The 
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solution  for  linear  and  angular  velocity  over  time  was  calculated  as  the  mass  was  reduced  at  a 
constant  rate.  At  time  t  =  0.5  the  mass  reduction  was  ceased.  The  results  in  figures  5-2  and 
5-3  show  linear  momentum,  P ,  and  angular  momentum,  H,  is  conserved  in  the  calculation. 
In  each  figure,  the  horizontal  line  is  the  constant  momentum  and  the  slanted  line  is  either  the 
mass  (figure  5-2)  or  the  moment  of  inertia  (figure  5-3).  The  parabolic  line  is  either  the  linear 
velocity  (figure  5-2)  or  the  angular  velocity  (figure  5-3). 


Figure  5-2:  Conservation  of  Linear  Momentum 

5.1.2  Validation  of  Compart ment  Flooded  Volume  and  Moment  of  Inertia 
Calculation 

The  flooded  compartment  module  added  to  LAMP  calculates  volume,  mass  center,  and  moment 
of  inertia.  The  module  also  calculates  a  final  hydrostatic  position  after  flooding  is  complete  in 
any  number  of  specified  compartments.  These  calculations  were  validated  by  running  them  with 
a  rectangular  hull  geometry  input.  Separate  calculations  were  made  with  the  same  rectangular 
hull  using  different  calculation  methods  and  software.  Results  from  the  LAMP  module  agreed 
with  results  from  the  independent  method  to  within  99%. 
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Figure  5-3:  Conservation  of  Angular  Momentum 

5.1.3  Validation  of  Flux  Difference  Splitting  Method 

A  bore  propagation  problem  is  solved  to  validate  the  flux  difference  splitting  method  used  for 
the  water  motion  on  deck  problem.  In  particular,  this  problem  validates  the  homogeneous 
form  of  the  difference  scheme,  equation  4.6.  The  initial  conditions  are  a  Riemann  Problem 

v  =  0.2667  ]  v  =  1.6 

>  y  <  0  and 

A  =  10.8  A  =  1.8 


where  the  velocity  units  are  tyi/s  and  the  water  height  units  are  m.  The  exact  solution  to  the 
problem  for  water  depth  has  been  given  by  Stoker  in  [30] : 

10.8  for  y  <  — 10.0< 

55^(20.84  -  |)2  for  -  lO.Of  <  y  <  0.45t 
4.716  for  0.45i  <  y  <  10.7 1 

1.8  for  y  >  10. 7t 

The  exact  solution  shows  that  a  bore  travels  in  the  positive  y-direction  with  a  velocity  of  10.7 
m/s,  and  a  rarefaction  wave  propagates  in  the  negative  y-direction  at  a  speed  of  10.0  m/s. 
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Figure  5-4  shows  the  numerical  results  (dashed  line)  and  the  exact  solution  (solid  line)  for  the 
water  depth.  The  numerical  results  match  well  with  the  exact  solution.  Because  the  numerical 
scheme  is  first-order,  it  has  some  difficulty  handling  the  shock  front.  This  is  evident  in  lost 
resolution  and  that  the  shock  front  travels  slightly  faster  than  the  exact  speed.  This  can  be 
corrected  with  a  higher  order  scheme. 


Figure  5-4:  Bore  Propagation  at  t  =  0.5  seconds 

A  second  problem  was  used  to  validate  the  flux  difference  splitting  method  when  the  ship  is 
undergoing  motion.  Mathematically  this  is  the  condition  where  the  right  hand  side  of  equation 
4.1  is  non-zero.  This  problem,  then,  is  a  test  on  the  numerical  scheme  used  for  equation  4.12. 
The  problem  is  one  solved  by  Huang  and  Hsiung  in  [12]  for  sloshing  inside  a  deck  well.  The 
computations  are  carried  out  for  the  deck  flow  excited  by  roll  motion.  The  deck  is  1  meter 
wide,  oscillating  about  a  pivot  which  is  0.522  meters  above  the  deck,  and  the  roll  amplitude  is 
0.067  rad.  The  mean  water  depth  is  6  cm.  The  primary  resonant  frequency  of  the  shallow- 
water  motion  inside  the  deck  well  is  u0  =  2.41  rad/  sec.  The  numerical  results  in  figure  5-5 
are  for  a  rolling  frequency  of  u>  —  0.9  rad j  sec,  well  below  resonance,  to  show  the  water  behaves 
as  a  ’’horizontal  surface.”  The  rolling  frequency  was  then  increased  to  twice  the  resonance 
frequency  and  results  plotted  in  figure  5-6.  At  this  frequency,  standing  waves  are  formed  with 
the  wave  length  approximately  equal  to  the  deck  width.  In  each  figure,  plots  were  made  of  the 
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free  surface  at  several  time  instances  to  illustrate  the  motion. 


Figure  5-5:  Shallow  Water  Sloshing  Below  Resonant  Frequency 


5.2  Results  for  Flooding  Models 

This  section  provides  some  results  from  the  compartment  flooding  model  that  was  added  to 
LAMP.  In  all  cases,  unless  noted  otherwise,  the  LAMP-2  formulation  was  used.  The  results 
include  effects  on  roll  and  pitch  motions  from  a  flooded  volume  and  examples  of  inaccuracies 
that  arise  due  to  a  shifting  eg  and  due  to  a  changing  mean  body  boundary  position  when  linear 
hydrodynamics  are  used.  The  ship  model  that  was  run  in  LAMP  for  these  results  is  a  CG47 
Class  Cruiser.  A  steady  forward  speed  of  10  knots  was  used  for  all  the  calculations.  This  ship 
model  was  used  on  the  recommendation  of  SAIC  due  to  its  excellent  past  performance. 

All  plots  of  ship  motion  were  made  in  non-dimensional  units.  If  the  ship  has  length  L  then 
linear  motions  are  made  non-dimensional  by  dividing  by  L  and  time  is  made  non-dimensional 
by  dividing  by  •  Angular  motions  are  plotted  in  radians. 
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Figure  5-6:  Shallow  Water  Sloshing,  Twice  Resonance  Frequency 

5.2.1  Roll  Motion  Results 

The  first  results,  figures  5-7  and  5-8,  show  the  effect  of  flooding  and  sloshing  on  roll  motion. 
Beam  seas  with  a  single  sinusoid  wave  component  were  used  to  induce  ship  roll  motions.  A 
large  amidships  compartment  with  length  equal  to  20%  of  ship  length  was  created  for  flooding. 
The  bulkhead  locations  were  at  0.4  and  0.6.  Figure  5-7  illustrates  roll  motion  of  the  intact  hull 
(solid  line)  and  of  the  flooded  ship  (dashed  fine).  For  the  flooded  ship,  the  problem  was  started 
with  the  ship  in  its  final  flooded  position  from  uncontrolled  flooding  into  the  compartment  and 
sloshing  was  not  included  in  the  calculation.  The  linear  hydrodynamics  and  sinusoid  seaway 
cause  the  motion  to  reach  steady  state  after  the  starting  transients  pass.  The  roll  motions  of 
the  flooded  ship  are  much  reduced  compared  to  the  intact  ship  due  to  the  flooded  mass  having 
lowered  the  ship’s  eg.  Sloshing  effects  were  then  added  to  the  calculation.  Figure  5-8  compares 
roll  motion  of  the  flooded  ship  without  sloshing  (solid  fine)  and  with  sloshing  (dashed  line). 
The  sloshing  causes  a  roll  moment,  the  effects  of  which  can  clearly  be  seen  in  the  increased  roll 
amplitude. 
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5.2.2  Vertical  Motion  Results 

The  next  set  of  results,  figures  5-9  and  5-10,  show  the  effect  of  flooding  on  vertical  motion  and 
relative  elevation  at  the  bow.  Head  seas  with  a  single  sinusoid  wave  component  were  used 
to  induce  pitch  and  heave  motions.  The  wavelength  was  set  at  1.5  times  the  ship  length  to 
excite  pitch  motion.  Wave  amplitude  was  17  feet  to  approximate  about  a  sea  state  6.  These 
extremely  harsh  conditions  were  intentional  to  create  a  large-  relative  elevation  at  the  bow  for 
water  shipping.  A  forward  compartment  with  bulkhead  locations  0.1  and  0.25  was  created  for 
flooding.  Figure  5-9  illustrates  pitch  motion  of  the  intact  hull  (solid  line)  and  of  the  flooded 
ship  (dashed  line).  For  the  flooded  ship,  the  problem  was  started  with  the  ship  in  its  final 
flooded  position  from  uncontrolled  flooding  into  the  compartment.  The  trim  of  the  flooded 
ship  is  indicated  in  the  initial  pitch  angle.  Sloshing  was  not  included  in  any  of  the  pitch 
calculations.  The  effects  from  the  longitudinal  moment  induced  by  the  flooded  water  and  its 
moment,  of  inertia  can  be  seen  in  the  flooded  ship  motions.  Notice  that  the  flooded  ship  spends 
a  greater  amount  of  time  than  the  intact  ship  with  the  bow  pitched  downwards  (positive  pitch 
angle).  In  such  a  condition  shipping  water  is  more  likely  to  occur  which  is  why  a  green  water 
model  is  also  necessary  for  LAMP  to  be  used  as  a  damage  stability  prediction  tool.  There 
is  little  difference  in  the  intact  and  flooded  ship  heave  motions  for  the  same  flooding  scenario. 
Figure  5-10  illustrates  this  motion  for  the  intact  ship  (solid  line)  and  flooded  ship  (dashed  line). 
The  motions  are  similar  but  the  flooded  ship’s  vertical  position  is  offset  due  to  its  lower  eg. 

The  relative  height  of  the  bow  deck  edge  above  the  instantaneous  free  surface  1  is  calculated 
for  water  shipping  and  can  be  plotted  for  the  pitch  motions  in  figure  5-9  This  relative  bow 
height  is  shown  in  figure  5-11  for  the  intact  ship  (solid  line)  and  the  flooded  ship  (dashed  line). 
Shipping  water  occurs  during  a  negative  relative  bow  height  (positive  relative  elevation).  The 
time  duration  and  magnitude  of  each  shipping  water  event  is  greater  for  the  flooded  ship. 

5.2.3  Loss  of  Accuracy  Examples 

Possible  causes  for  loss  of  accuracy  in  the  LAMP  flooding  simulation  calculation  have  been 
discussed  and  examples  are  now  provided.  The  next  set  of  results  shows  that  a  changing  mean 

xThe  term  "relative  elevation”  was  used  in  the  green  water  model  discussion  as  meaning  the  free  surface 
height  minus  the  deck  height,  the  relative  bow  height  is  just  the  negative  value  of  the  relative  elevation. 
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Figure  5-8:  Sloshing  Effects  on  Roll  Motion 
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Figure  5-9:  Pitch  Motion  Intact  and  Flooded  Ships 


Figure  5-10:  Heave  Motion  Intact  and  Flooded  Ship 
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body  boundary  position  affects  linear  calculation  accuracy  and  also  that  a  shifting  eg  introduces 
errors  in  the  angular  motion  calculation.  Figure  5-12  plots  pitch  motion  to  illustrate  loss  of 
linear  calculation  accuracy.  For  each  curve  in  the  figure  the  problem  was  started  with  the  ship 
in  its  final  flooded  position  from  uncontrolled  flooding  in  the  forward  compartment  bounded 
by  the  0.1  and  0.25  bulkheads.  The  seaway  conditions  were  the  same  as  the  conditions  for  the 
figure  5-9  simulation.  The  solid  curve  is  the  accurate  linear  hydrodynamic  calculation;  it  used 
the  final  flooded  hydrostatic  position  as  the  mean  body  boundary  position.  The  dashed  curve 
used  the  ship’s  intact  (non-flooded)  mean  body  position  for  the  linear  hydrodynamic  calculation. 
This  example  shows  the  calculation  with  the  incorrect  mean  body  position  overestimates  the 
peak  positive  pitch  by  approximately  0.5  degrees. 

Figure  5-13  illustrates  the  effect  on  angular  motion  calculation  due  to  a  shift  in  eg.  For 
each  plot  the  problem  was  started  with  the  ship  in  its  final  flooded  position  from  uncontrolled 
flooding  in  the  forward  compartment  bounded  by  the  0.1  and  0.25  bulkheads.  This  flooded 
condition  moves  eg  forward  by  3.8%  and  down  0.23%  of  ship  length  from  its  initial  intact 
position.  The  seaway  conditions  were  the  same  as  the  conditions  for  the  figure  5-9  simulation. 
The  solid  curve  is  the  pitch  motion  where  the  ship’s  eg  was  moved  prior  to  the  start  of  the 
calculation  to  account  for  the  flooded  water.  The  dashed  curve  is  the  same  calculation  except 
that  the  ship’s  eg  was  only  moved  longitudinally  to  account  for  the  flooded  water.  It  was  not 
moved  vertically. 

5.2.4  Progressive  Flooding  Results 

The  last  flooding  simulation  result  is  an  example  of  progressive  flooding.  The  seaway  conditions 
were  the  same  as  the  conditions  for  the  figure  5-9  except  that  the  wave  amplitude  was  reduced 
by  two  thirds.  The  problem  was  started  with  the  ship  in  its  final  flooded  position  from 
uncontrolled  flooding  in  the  forward  compartment  bounded  by  the  0.1  and  0.25  bulkheads.  At 
time  equal  to  12,  indicated  by  ’’start”  in  figure  5-14,  flooding  at  a  constant  rate  was  initiated 
into  a  compartment  bounded  by  bulkheads  0.25  and  0.4.  At  time  equal  to  24,  indicated  by 
’’stop”  in  figure  5-14,  flooding  ceased.  Figure  5-14  plots  the  pitch  and  heave  motions,  the 
heave  motion  is  the  lower  curve  in  the  figure.  Figure  5-15  plots  the  relative  bow  height  during 
the  progressive  flooding  simulation.  Note  that  shipping  water  events  do  not  occur  until  the 
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progressive  flooding  commences. 

As  a  comparison,  this  progressive  flooding  simulation  was  run  a  second  time  using  the  fully 
non-linear  LAMP-4  formulation.  The  linear  and  non-linear  calculation  results  for  pitch  motion 
and  relative  bow  height  are  plotted  in  figures  5-16  and  5-17.  The  non-linear  calculation  results 
are  plotted  as  solid  fines,  the  linear  as  dashed  fines.  The  results  for  the  two  calculation  methods 
have  close  agreement.  As  before,  however,  the  calculation  with  the  less  accurate  mean  body 
boundary  position  (linear  calculation)  overestimates  the  pitch  motion. 

5.3  Results  for  Green  Water 

This  section  provides  some  results  from  the  green  water  model  that  was  added  to  LAMP.  In 
all  cases  the  LAMP-2  formulation  was  used.  The  results  include  effects  on  vertical  motions 
from  green  water  and  plots  of  green  water  mass  and  local  deck  loads.  A  CG47  Class  Cruiser  at 
a  steady  forward  speed  of  10  knots  was  used  for  the  calculations.  Plots  of  ship  motion,  mass, 
and  mass  centers  are  made  in  non-dimensional  units.  Because  the  CG47  is  such  a  massive  ship, 
the  density  of  the  green  water  used  in  calculations  for  all  of  the  results  below  was  increased  by 
a  factor  of  five.  This  was  done  to  make  the  green  water  mass  and  its  effects  on  ship  motion 
stand  out. 

Head  seas  with  a  single  sinusoid  wave  component  were  used  .  The  wavelength  was  set  at 
1.5  times  the  ship  length  to  excite  pitch  motion.  Wave  amplitude  was  14  feet.  This  amplitude 
was  determined  through  trial  and  error  to  provide  about  a  10  foot  relative  elevation  for  the 
water  shipping  problem.  Relative  elevations  higher  than  this  required  increased  calculation 
time  because  smaller  time  steps  were  needed  to  solve  the  shallow  water  problem.  Also,  for 
relative  elevations  much  larger  than  10  feet  it  was  found  that  large  time-dependent  mass  and 
moment,  of  inertia  terms  were  calculated  that  caused  the  dynamic  equations  of  motion  solution 
to  break  down. 

Figure  5-18  shows  the  effect  of  green  water  on  pitch  motion  This  figure  plots  the  relative 
bow  height  and  pitch  motion  of  a  LAMP  calculation  with  no  green  water  effects  (solid  lines) 
and  with  green  water  effects  (dashed  fines).  The  relative  bow  height  is  the  curve  with  smaller 
negative  amplitudes.  The  increased  angular  inertia  of  the  ship  with  the  green  water  delays  its 
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Figure  5-12:  Effects  of  Linear  Hydrodynamics  and  Flooding 
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Figure  5-13:  Effect  of  Center  of  Gravity  Shift  on  Pitch  Motion 


Figure  5-14:  Progressive  Flooding  Pitch  and  Heave  Motion 
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motion  response  to  the  incoming  waves  and  slightly  increases  the  maximum  amplitude  of  the 
positive  pitch  motion.  Figure  5-19  shows  the  effect  of  green  water  on  heave  motion.  Heave 
motion  is  slightly  increased  in  the  motion  calculation  that  includes  green  water  (dashed  line). 
Again,  the  green  water  density  was  artificially  increased  by  a  factor  of  five  for  these  results. 
For  smaller  ships  and  craft  the  green  water  would  have  a  much  larger  effect  on  pitch  motion. 
Green  water  simulations  were  not  run  for  beam  seas  but  roll  amplitudes  would  probably  also 
be  enhanced  due  to  a  greater  tendency  of  a  ship  to  roll  than  pitch. 

The  calculation  for  V  T  in  the  water  shipping  problem  was  performed  using  Methods  I,  II, 
and  III.  Figure  5-20  plots  relative  bow  height  (the  sinusoidal  curve)  and  green  water  mass  on 
deck  for  two  water  shipping  events  using  Method  I  for  water  shipping  problem.  The  green  water 
mass  was  normalized  by  dividing  it  by  the  ship  mass.  A  second  calculation  using  Method  II 
for  the  water  shipping  problem  showed  no  discernible  differences  from  the  Method  I  calculation 
results.  The  amount  of  green  water  mass  shipped  on  board  in  each  calculation  was  the  same. 

The  Method  III  calculation  for  shipping  was  based  on  conserving  momentum  of  the  incoming 
wave.  Figure  5-21  plots  the  green  water  mass  on  deck  using  Method  III  (dashed  line)  and  the 
green  water  mass  using  Method  I  (solid  line).  The  Method  IH  calculation  for  shipping  water 
velocity  should  be  the  method  normally  used  because  it  results  in  the  most  green  water  mass 
on  deck;  it  provides  the  most  conservative  estimate  of  water  shipping. 

The  green  water  model  solved  the  water-motion-on-deck  problem  in  strips  of  width  dely 
assuming  a  two-dimensional  free  surface.  This  approach  prevents  transverse  water  flow  so  that 
no  green  water  is  able  to  fall  off  the  port  and  starboard  sides  of  the  weatherdeck.  Therefore 
the  time-duration  that  the  green  water  mass  is  calculated  to  be  on  the  weatherdeck  is  longer 
than  the  time-duration  if  a  three-dimensional  free  surface  were  solved.  For  three-dimensional 
calculations,  the  green  water  mass  curve  would  have  about  the  same  peak  value  but  would  go 
to  zero  after  each  water  shipping  event  in  a  shorter  amount  of  time.  It  would  have  more  of  a 
spiked  appearance. 

The  green  water  mass  center  location  and  local  deck  loads  are  also  presented.  Figure  5- 
22  plots  mass  and  mass  center  for  two  water  shipping  events.  The  weatherdeck  longitudinal 
coordinate  between  0.4  and  0.5  corresponds  to  the  bow  area.  The  aft  end  of  the  weatherdeck 
has  coordinate  0.25.  Figure  5-23  plots  local  deck  pressure  due  to  green  water  and  pitch  angle 
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(sinusoidal  plot).  The  dashed  line  is  pressure  for  a  location  at  the  bow,  the  solid  line  is  pressure 
for  a  location  aft  of  the  bow  a  distance  0.06  (non-dimensional).  The  green  water  pressures  are 
normalized  by  dividing  them  by  the  hydrostatic  head  of  the  green  water,  pg(,  where  £  is  the 
green  water  depth.  This  figure  shows  that  the  green  water  deck  pressure  increases  by  almost 
a  factor  of  two  over  hydrostatic  during  peak  accelerations  in  pitch  motion. 
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Figure  5-18:  Effects  of  Green  Water  on  Pitch  Motion 
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Figure  5-19:  Effects  of  Green  Water  on  Heave  Motion 


Figure  5-20:  Relative  Bow  Height  and  Green  Water  Mass  on  Deck 
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Figure  5-21:  Mass  on  Deck  Using  Different  Shipping  Water  Velocities 


Figure  5-22:  Green  Water  Mass  on  Deck  and  Mass  Center 


TO 
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Chapter  6 


Conclusions 


6.1  Discussion  and  Recommendations 

This  research  investigated  the  addition  of  models  for  compartment  flooding  and  green  water  to 
the  LAMP  time-domain  ship  motion  program  so  that  it  could  be  used  as  a  damaged  stability 
prediction  tool.  The  approach  was  to  treat  both  events,  the  flooded  water  and  the  green 
water  on  the  weatherdeck,  as  a  change  in  the  ship  rigid  body  mass.  Time-dependent  mass 
and  mass  moment  of  inertia  terms  were  calculated  for  the  water  and  incorporated  into  the 
dynamics  equations  of  motion.  The  results  from  calculations  with  the  models  showed  that  ship 
motions  are  indeed  affected  by  the  flooding  and  green  water  events  and  that  these  motions  can 
be  estimated  with  the  time-domain  program.  It  was  also  found  that  for  a  large  ship  green 
water  has  little  effect  on  ship  motions  but  that  it  can  cause  a  significant  increase  in  local  loads 
on  the  weatherdeck  and  superstructure.  In  short,  the  models  developed  by  this  thesis  along 
with  the  LAMP  program  can  be  used  as  a  damaged  stability  prediction  tool. 

A  major  problem  with  running  the  LAMP  3-D  nonlinear  formulation  (LAMP-4)  is  the 
extensive  computation  time,  especially  to  run  a  full  flooding  simulation  over  a  long  timeline. 
Unless  ship  motions  are  very  large,  however,  the  nonlinear  calculation  is  probably  unnecessary. 
A  comparison  of  calculation  results  using  the  nonlinear  LAMP-4  formulation  and  the  linear 
LAMP-2  formulation  shows  that  the  LAMP-2  results  provide  adequate  accuracy  in  the  ship 
motion  calculation  for  estimating  damaged  stability.  For  best  results  with  the  linear  calculation 
for  a  flooding  simulation,  an  appropriate  mean  body  boundary  position  needs  to  be  selected 


72 


by  careful  consideration  of  the  ship’s  initial  and  final  flooded  position.  Experience  with  the 
hydrodynamic  properties  of  the  hull  in  question  would  help  in  determining  what  body  position 
to  linearize  about. 

Several  methods  of  calculating  the  relative  velocity  for  the  water  shipping  problem  were 
attempted.  It  was  concluded  that  the  best  method  was  Method  III  which  conserved  the,  wave 
particle  momentum. 

6.2  Problems  Encountered 

6.2.1  Computational  Difficulties  During  Rapid  Changes  in  Mass  and  Mass 
Distribution 

Time-dependent  mass  and  mass  moment  of  inertia  was  included  in  equation  2.45  for  the 
dynamics  equations  of  motion.  A  backward  difference  scheme  was  used  to  calculate  the  time- 
derivative  terms.  Say  fi(£)  stands  for  either  the  mass  matrix  or  the  mass  moment  of  inertia 
tensor,  then  =  Bih 1) At)  represents  the  time-derivative  for  the  ith,jth 

component  of  the  matrix  or  tensor.  During  simulations  where  these  derivative  terms  were 
large,  for  example  water  shipping  with  relative  elevations  greater  than  about  10  feet  or  a  large 
flooded  volume  undergoing  sloshing  with  roll  amplitudes  greater  than  25  degrees,  it  was  found 
the  dynamic  solution  method  became  unstable.  Changing  the  time  derivative  to  ^  = 

n(t ,j,n At) ,(n-2)At)  provided  more  of  an  averaged  estimate  for  the  derivative  and  sometimes 
kept  the  dynamic  solution  method  stable.  Still,  some  calculations  with  large  mass  and  mass 
moment  of  inertia  time-derivative  values  of  were  unable  to  be  completed  due  to  instability. 
More  elaborate  methods  of  estimating  the  derivatives  were  not  investigated. 

6.2.2  Selection  of  the  Time  Discretization  for  the  Flux  Difference  Splitting 
method 

For  a  given  discretization,  the  stability  of  the  flux  difference  splitting  method  used  to  solve  the 
water  motion  on  deck  problem  is  a  function  of  the  water  depth  and  particle  velocity.  As  either 
of  these  quantities  becomes  too  large  the  scheme  can  become  unstable  unless  the  length  of  the 
time  step  is  decreased.  Smaller  time  steps,  however,  slow  the  LAMP  calculation  and  make  it 
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computationally  expensive  to  get  through  a  complete  damaged,  scenario  or  flooding  timeline. 
Also,  the  maximum  relative  elevation  for  water  shipping,  which  governs  the  green  water  depth, 
may  not  be  known  prior  to  completing  the  calculation.  If  in  this  case  the  maximum  relative 
elevation  is  underestimated  during  calculation  setup  then  the  calculation  may  be  unstable;  if 
the  elevation  is  overestimated  during  setup  then  the  time  step  will  be  smaller  than  necessary. 

For  a  large  ship  where  green  water  has  little  effect  on  ship  motions,  the  LAMP  calculation 
could  be  uncoupled  from  the  green  water  calculation  by  running  the  green  water  calculation  as 
a  post  process  using  data  generated  by  the  main  LAMP  calculation.  The  LAMP  generated 
data  could  be  analyzed  in  order  to  optimally  setup  the  green  water  post  process  calculation. 
Of  course,  for  a  small  ship  in  large  seas  green  water  will  affect  motions  and  the  green  water 
calculation  should  remain  coupled  with  the  main  LAMP  calculation. 

6.2.3  Calculating  Relative  Velocity  for  the  Water  Shipping  Problem 

Accurate  calculation  of  relative  velocity  for  the  water  shipping  problem  requires  more  research. 
Sea  spectrum  are  typically  described  using  a  summation  of  gravity  waves  each  of  which  solves 
the  linearized  free  surface  boundary  condition  and  with  the  wave  particle  velocities  specified 
on  the  undisturbed  plane  of  the  free  surface  and  below.  To  accurately  calculate  the  relative 
velocity  when  shipping  water,  the  wave  particle  velocity  for  the  portion  of  the  wave  that  is 
instantaneously  higher  than  the  ship  bulwark  or  deck  edge  must  be  known. 

If  the  water  shipping  event  occurs  when  the  deck  edge  or  bulwark  is  above  the  undisturbed 
free  surface,  or  with  large  amplitude  waves,  then  linear  wave  theory  may  not  provide  an  accurate 
value  for  water  particle  velocity.  Shipping  water  is  inherently  a  nonlinear  problem  but  has  been 
modeled  in  this  thesis  with  linear  assumptions.  Stokes  expansions,  which  are  nonlinear  solutions 
for  plane  waves  based  on  systematic  power  series  in  the  wave  amplitude,  may  be  a  method  of 
describing  the  seaway  to  introduce  nonlinear  water  particle  velocity  into  the  water  shipping 
calculation. 

Finally,  hydraulic  flow  was  used  in  the  relative  velocity  calculation  but  this  also  an  approx¬ 
imation  of  the  actual  velocities  that  occur  at  the  wave  and  deck  edge  interface.  Some  water 
shipping  calculations  were  performed  by  omitting  the  hydraulic  velocity  term  but  the  results 
did  not  look  very  physical.  Too  little  green  water  mass  entered  onto  the  weatherdeck.  It  is 
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concluded  that  in  addition  to  the  terms  for  particle  and  deck  edge  velocity  a  third  term  such 
as  the  hydraulic  term  is  need  for  accurate  relative  velocity  calculation.  However,  the  hydraulic 
method  may  not  be  the  correct  way  to  calculate  the  term.  Use  of  wave  phase  or  group  velocity, 
or  applying  dam-breaking  results  from  [29],  may  provide  an  alternate  method  of  calculating 
the  third  term  to  give  a  more  accurate  water  shipping  model.  Theoretical  uncertainties  in 
calculating  the  relative  water  velocity  explain  why  the  constant  C0  must  be  used  in  equation 
4.30  to  adjust  theoretical  estimates  to  match  experimental  results. 

6.3  Recommendations  for  Future  Research 

Although  this  research  has  demonstrated  including  models  for  flooding  and  green  water  into 
LAMP  so  that  it  can  be  used  as  a  damaged  stability  prediction  tool,  much  more  work  should 
be  done  to  develop  the  models.  Recommendations  include: 

•  The  green  water  model  currently  consists  of  longitudinal  strips  on  the  weatherdeck  where 
two-dimensional  free  surface  calculations  are  performed  to  solve  the  water  flow  on  deck 
problem.  This  method  prevents  transverse  flow  of  water  which  does  not  allow  water  to  fall 
off  the  sides  of  the  weatherdeck.  The  result  is  that  the  calculated  green  water  mass  on  the 
ship  becomes  artificial  and  remains  on  the  ship  longer  than  it  should.  A  three-dimensional 
free  surface  green  water  model  should  be  developed  to  correct  this  condition. 

•  Some  experimentation  has  been  performed  on  water  shipping  such  as  that  reported  by 
Grochowalski  in  [11].  These  experiments  should  be  repeated  in  a  computer  simulation 
with  the  LAMP  models  in  order  to  validate  the  models.  This  work  would  also  help 
in  determining  a  proper  value  for  the  constant  in  the  water  shipping  mass  introduction 
equation  4.30. 

•  A  rich  area  for  future  work  is  to  examine  the  nonlinear  aspects  of  the  water  shipping 
problem  and  formulate  a  more  theoretical  approach  for  calculating  relative  velocity  in  the 
water  shipping  problem.  The  scope  of  this  work  should  also  look  at  calculation  of  the 
free  surface  elevations  when  considering  the  pile-up  that  occurs  from  slamming. 
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•  Damaged  model  testing  was  reported  in  references  [14],  [15],  and  [16].  Work  should  be 
performed  to  see  if  these  results  can  be  repeated  in  LAMP  simulations. 

•  Finally,  structural  loads  on  the  ship  while  in  a  damaged  condition  have  not  been  addressed 
in  this  thesis  except  for  illustrating  that  the  local  deck  loads  will  increase  in  the  green 
water  problem.  The  flooding  and  green  water  models  should  be  expanded  to  include 
effects  on  main  girder  and  local  structural  loads. 
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Appendix  A 


Moment  of  Inertia  Tensor 
Calculations 


This  appendix  provides  notes  on  the  moment  of  inertia  tensor  calculations  for  a  flooded  com¬ 
partment  volume.  Figure  A-l  shows  a  parallelepiped  with  sides  of  length  a,  6,  and  c  as 
indicated. 


Figure  A-l:  Parallelepiped 


The  volume  moment  of  inertia  of  the  parallelepiped  about  its  center  is }IXx 


=  ^(b2+c 2), 
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lyy  =  «|(a2  +  ^  and  izz  =  ^(a2  +  b2).  The  moment  of  inertia  tensor,/,  is  written  as, 

4*  0  0 
0  lyy  0 

0  0  Izz 

The  parallel  axis  theorem  can  be  used  to  refer  7  to  the  primed  coordinate  system  in  figure 
A-l.  The  coordinates  ( A ,  B,  C).  are  the  center  of  the  parallelepiped  in  the  primed  frame.  The 
equation  for  the  parallel  axis  theorem  is 

B2  +  C 2  -AB  -AC 

7/  =  7  +  abc  -AB  A2  +  C 2  -BC 

-AC  -BC  A2  +  B2 

If  the  parallelepiped  moment  of  inertia,  J,  was  obtained  about  a  point  other  than  its  center, 
then  there  would  be  additional  terms  in  the  parallel  axis  equation. 

The  rotational  transformation  of  the  moment  of  inertia  tensor  refers  the  tensor  to  a  co¬ 
ordinate  system  that  is  related  to  the  first  through  a  rotational  transformation  matrix,  C. 
This  matrix  performs  the  sane  function  as  the  euler  angle  transformation  matrix,  L,  discussed 
in  Chapter  2.  Say  the  moment  of  inertia  tensor  of  a  body  I  is  referred  to  the  unprimed 
coordinate  system  in  figure  A-2  and  C  is  the  rotational  transformation  matrix  between  the 
twocoordinate  systems.  The  transformation  of  I  is  performed  by, 

d  =  cTct 
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Appendix  B 


Details  of  Riemann  Problem 
Calculations 


In  reference  (stoker),  Stoker  demonstrates  how  the  Riemann  problem  may  be  solved  when  the 
fluid  is  initially  at  rest  on  the  right-hand  side  of  the  dam.  In  order  to  use  some  of  the  Stoker 
results  in  the  below  calculations,  an  axis  system  will  be  adopted  that  moves  at  constant  velocity. 
Use  of  the  Stoker  results  is  justified  by  showing  that  the  solution  to  the  dam  breaking  problem 
is  invariant  with  respect  to  the  moving  axis  system. 

Moving  Axis  System 

Figure  B-l  illustrates  initial  conditions  for  the  Riemann  problem  in  stationary  frame  y  while 
figure  B-2  illustrates  the  same  problem  in  reference  frame  y'  translating  at  steady  velocity,  Vo- 
The  two  reference  frames  are  related  as  follows,  where  c,  the  wave  propagation  speed,  equals 
y/gT. 


y  =  y'  +  V0t 

(B-l) 

V/(yf,t)  =  V(yf  +  V0t,t)-V0 

(B-2) 

A/(y/,  t)  =  A  (yf  +  V0t,  t) 

(B.3) 
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Figure  B-l:  Initial  Conditions  for  the  Riemann  Problem 
cf(yf,  t)  =  c(yt  +  V0t,  t ) 


(BA) 


From  reference  (stoker)  the  shallow  water  wave  equations  can  be  formulated  in  terms  of  the 
velocity,  V,  and  propagation  speed,  c.  These  equations  for  the  stationary  frame  are 


(B.5) 

(B.6) 

which  can  be  added  or  subtracted  from  each  other  to  form  characteristic  equations 


{l+(F+c)|}(l,+2c>=0  (BJ) 

{I+(v'c)|;}(l'_2c)=0  (B-8) 

Equations  B.7  and  B.8  state  that  the  function  (V  +  2c)  is  constant  for  a  point  moving 
through  the  fluid  with  the  velocity  (V  +  c)  and  that  the  function  (V  —  2c)  is  constant  for  a 
point  moving  through  the  fluid  with  the  velocity  (V  —  c).  There  are  two  sets  of  curves,  Ciand 
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Figure  B-2:  Initial  Conditions  for  Riemann  Problem  in  Steady  Velocity  Frame 


C2  called  characteristics,  which  are  the  solution  curves  of  the  ordinary  differential  equations 

Ci  :  ^  =  V  +  c  (B.9) 

Cut 

Cl  :  %  =  V-C  (B.10) 

at 


Along  characteristic  C\,  V  +  2c  is  constant  and  along  C2,  V  +  2c  is  constant.  Similarly,  the 
shallow  water  wave  equations  for  the  translating  frame  are 


dVt 

dt 


+  Vf 


dvt 

dyf 


+  2  d 


dd 

dyf 


=  0 


(B.ll) 


dd  n.r  dd  ,  dVf  . 
2—  +  2 Vi-t-—  +  —  9 

dt  dyf  dyf 


(B.12) 


substituting  equations  B.2  and  B.3  into  equations 


v°8£+w+<-v-v°^+2% 


=  0 


y=y/+Vot 


(B-13) 


=0 


y=yf+Vot 


(B.14) 
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Equations  B.13  and  B.14  can  be  simplified  to 


dV  „dV 
dt  +V  dy 


n  dc 
+  2  c— 
dy 


y=yi+Vot 


y=y'+Vot 


0 

0 


(B.15) 

(B.16) 


which  are  the  same  equations  as  B.5  and  B.6.  Thus  the  translating  frame  has  the  same 
solution  curves,  characteristics,  as  the  stationary  frame  which  shows  the  solution  to  B.5  and 
B.6  are  invariant  with  respect  to  axes  moving  with  constant  velocity. 


Relations  From  Stoker 

The  following  results  from  reference  (stoker)  will  be  used  in  the  details  of  the  Riemann  Problem 
calculation.  Referring 

to  Figure  B-3,  for  a  bore  advancing  at  velocity  £,  the  equations  for  mass  and  momentum 
conservation  are 


1 

g 

o 

II 

<77 

1 

£ 

rH 

-3. 

(B.17) 

p\\V\  (Vi  —  £)  —  pAoVo  (Vo  —  £)  —  ~  ^i) 

(B.18) 
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equation  B.17  can  be  solved  for  the  shock  speed,  £ 


AiVj  -  A0V0 

Ai  —  Aq 


(B.19) 


equations  B.17  and  B.18  can  be  combined  to  solve  for  V\ 


V\  = 


g  (Aq  +  Ai)  (Ai  —  Aq)2 

2  AqAi 


+  V0 


(B.20) 


When  Vo  =  0,  the  following  relations  can  be  derived  from  equations  B.17  and  B.18 


Vi 

Co 


1 

Co 


(B.21) 


Calculations  for  Case  I 


(B.22) 


When  V2  +  2c2  >  Vo  +  2co,  the  values  for  the  C\  characteristics  of  the  high  side  of  the  dam  and 
the  low  side  of  the  dam  are  discontinuous  which  forms  a  shock,  represented  by  £/  in  figure  B-4. 
Velocity  and  water  height  values  are  primed  in  figure  B-4  due  to  the  coordinate  translating  at 
steady  velocity  Vo- 

Curved  characteristic  Ch  indicated  in  figure  B-4,  is  the  solution  curve  where  V3/  +  2 c3r  is 
constant.  Since  C\  intersects  dashed  line  I  in  zone  1,  then  C\  =  Vi f  +  2ci/  .  Therefore 
V3/  +  2c3/  =  Vi/-|-2ci/.  Also,  since  C\  intersects  dashed  line  II  in  zone  2,  then  C\  =  V2f  +  2c2>. 
Therefore  V3f  +  2c3t  =  V2t  +  2 c2f.  Then 


Vi/  +  2cx/  =  V2f  +  2c2/ 


(B.23) 


substituting  equation  B.21  for  Vj/  and  equation  B.22  for  C]/,  equation  B.23  can  be  written 
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Figure  B-4:  Case  I,  Coordinate  System  Translating  at  Vo 


next  substitute  R  —  ^  into  equation  B.24 


V2  /  +  2c2/ 
Co' 


(B.25) 


After  solving  equation  B.25  for  R,  the  solution  for  zone  1  can  be  determined  for  the  fixed 
axis  system  by  substituting  £/  =  £  —  Vq,  and  V\'  =  V\  —Vq. 


There  are  a  fan  pattern  of  straight  line  characteristics,  not  shown  in  figure  B-4,  that  fan 
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between  characteristic  I  and  characteristic  II  and  pass  through  characteristic  C\.  On  C 1 


Vst  +  2c3/  =  V2t  +  2c2t  (B.26) 

Also,  on  the  fan  of  straight  characteristics  in  zone  3 

^  =  £  =  VZ!  -  (B.27) 

at  t 

combining  equations  B.26  and  B.27  to  eliminate  c3/  gives 

(a28) 

finally  determine  the  zone  3  velocity  solution  in  the  fixed  axis  system  by  substituting  V2f  = 
V2  -  Vo,  V3f  =  V3  -  V0,  and  f  =  f  -  V0. 

v*=i(i  -v«+^-2^+^)+v» 

Similar  to  the  solution  for  V3.  c3  can  be  solved  for  by  combining 
to  eliminate  V3t  giving 

c3/  =  -  (v2!  +  2c2l  — 
substituting  A 3/  =  (C3J)  ■  into  B.29  gives 

U  =  +  f)2 

finally  determine  the  zone  3  height  solution  in  the  fixed  axis  system  by  substituting  V2r  = 
V2  —  Vo,  and  =  3L  _  Vo. 

*=ra{Vl+202- if 

Zone  3  is  bounded  by  characteristic  I  and  characteristic  II.  For  characteristic  I,  = 
W  -  ci/,  and  for  characteristic  II,  =  Vy  -  c2/.  Therefore  the  zone  3  solution  is  bounded 


equations  B.26  and  B.27 
(B.29) 

(B.30) 
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by 

(V2f  —  c2f)  t  <  yt  <  {V\i  —  ci/)  t  (B.31) 

or  converting  to  the  fixed  axis  system,  the  zone  3  solution  is  bounded  by 

{V2-c2)t<y<(Vi-ci)t  (B.32) 

The  above  solution  for  the  zone  3  velocity  are  for  a  velocity  that  increases  with  increasing  y. 
Thus,  the  Case  I  results  are  only  good  when  V2  <  V\  to  allow  for  an  increasing.^  If  V2  >  Vi, 
a  velocity  step-down  between  zones  2  and  1  would  form  a  shock  on  the  left  side  of  the  dam. 
Therefore,  after  solving  for  V\  make  the  following  check  to  see  if  V2  <Vi.  Equation  B.20  is 

substituted  for  V\  1 

y2  <  rg(Ao  +  Ai)(Ax  -AqH  2  +  (B.33) 

2  AoAi 

If  equation  B.33  is  not  met,  then  shocks  form  on  the  left  and  right-had  sides  of  the  dam 
and  the  Riemann  problem  solution  is  of  Case  II,  as  indicated  in  figure  B-5. 

Calculations  for  Case  II 


Figure  B-5:  Case  II,  Fixed  Coordinate  System 
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To  solve  Case  II,  equation  B.20  must  be  used  for  each  shock,  and  £2. 


Vi  = 


g  (Aq  +  Ai)  (Ai  —  Ap)2 
2  A0A1 


+  M, 


(B.34) 


V2  =  [£(Ai  +  A2)(A2-AiH  2  +  ^  (B.35) 

2  A1A2 

substituting  equation  B.34  into  equation  B.35  for  Vx  results  in  an  equation  that  can  be 
solved  for  Ai 

rr  rr  [ 5  (Ao  +  Ai )  ( Ai  -  Ao) 2 ]  2  g  (Ai  +  A2)  (A2  -  Ai)2 

" Fo "  [2  A^  J  [2  AXA2 

once  Aiis  determined  then  equation  B.34  can  be  used  to  calculate  V\.  The  two  shock  speeds 
are  found  using  equation  B.19 


A1V1  -  AqVq 
Ai  —  Aq 


and 


_  AaVa-A^i 

A2  -  Ai 


(B.37) 


Calculations  for  Case  III 

The  entering  argument  for  case  I  was  that  V2  +  2c2  >  Vo  +  2co-  When  this  condition  is  not 
met,  rarefaction  waves  instead  of  shock  waves  occur.  For  case  III,  V2  +  2c2  <  Vb  +  2co  which 
can  be  re-written  as  Vq  —  V2  >  2  (c2  —  Co)  or 


Vo  —  V2  >  2|co  —  c2| 


(B.38) 


The  solution  to  case  III  can  be  found  by  analyzing  characteristics  I  through  IV,  C\,  and 
C2  as  illustrated  in  figure  B-6. 

On  characteristic  C\,  Vz  +  2cz  is  constant.  Since  this  characteristic  intersects  zones  1  and  2 

V2  +  2c2  =  V1  +  2c1  (B.39) 
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Figure  B-6:  Case  III,  Fixed  Coordinate  Sysem 
On  characteristic  C%,  V4  -  2 c4  is  constant.  Since  this  characteristic  intersects  zones  1  and  0 

V0-2c0  =  V1-2c1  (B.40) 


adding  equations  B.39  and  B.40  gives 


Vi  =  — - - 1-  C2  —  co 


(B.41) 


subtracting  equations  B.39  and  B.40  gives 

C2  +  co  ,  V2  —  Vo 


ci  = 


+  • 


(B.42) 


since  Ai  =  equation  B.42  can  be  written  as 


C2  +  Co  V2  ~Vq 


(B.43) 


Zone  3  is  solved  similarly  to  the  zone  3  solution  in  case  I.  There  are  a  fan  pattern  of 
straight  line  characteristics,  not  shown  in  figure  B-6,  that  fan  between  characteristic  III  and 
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characteristic  IV  and  pass  through  characteristic  C\ .  On  C\ 

V3  +  2c3  =  V2  +  2c2  (B.44) 

Also,  on  the  fan  of  straight  characteristics  in  zone  3 

%  =  Z  =  V,-«  (B.45) 

dt  t 

combining  equations  B.44  and  B.45  to  eliminate  c3  gives 

(R46) 

Similarly,  solve  for  c3  by  combining  equations  B.44  and  B.45  to  eliminate  V3  giving 

<a  =  i  (vt  +  2c,  -  2 )  (B.47) 

substituting  A3  =  into  B.47  gives 

As  =  i(v2  +  2c2-f)2  (B.48) 

Zone  3  is  bounded  by  characteristics  III  and  IV.  For  characteristic  III,  %  —  V\—  cx,  and 
for  characteristic  IV,  ^  =  V2  —  02-  Therefore  the  zone  3  solution  is  bounded  by 

(V2-c2)t<y<(y1-c1)t  (B.49) 

For  Zone  4  there  are  a  fan  pattern  of  straight  line  characteristics,  not  shown  in  figure  B-6,  that 
fan  between  characteristics  I  and  II.  All  pass  through  characteristic  C2.  On  C2 

Vi-2ci  =  V0-2c0  (B.50) 

Also,  on  the  fan  of  straight  characteristics  in  zone  4 

^  =  f  =  y1+C4  (B.51) 

dt  t 
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combining  equations  B.50  and  B.51  to  eliminate  c4  gives 

(R52) 

Similarly,  solve  for  C3  by  combining  equations  B.50  and  B.51  to  eliminate  V4  giving 

Cl  =  i(f-V„+2 (B.83) 
substituting  A4  =  into  equation  B.53  gives 

Al=^d_Fo+2co)2 

Zone  4  is  bounded  by  characteristics  /  and  II.  For  characteristic  I,  -^  = 
characteristic  II,  ^  =  V\  +  cj.  Therefore  the  zone  4  solution  is  bounded  by 

(Vi  +c{)t  <  y  <  (Vo  +  co)t  (B.55) 

Calculations  for  Case  IV 

If  A0  =  0  then  the  solution  consists  of  solving  the  single  rarefaction  illustrated  in  figure  B-7. 
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'  A2 

V3\  1  NSl 
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Figure  B-7:  Case  IV,  Fixed  Coordinate  System 

The  solution  to  case  IV  can  be  found  by  analyzing  characteristics  I,  II,  and  C\.  There 
is  a  fan  pattern  of  straight  line  characteristics,  not  shown  in  figure  B-7,  that  fan  between 
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characteristics  I  and  II  and  pass  through  characteristic  G\.  On  C\ 


V3+2c3  =  V2  +  2c2  (B.56) 

On  characteristic  I,  the  depth  at  the  water  boundary  is  zero,  c3  =  0,  so  that  V3  =  V2  +  2 c2. 
This  means  that  the  boundaiy  of  the  water  moves  to  the  right  on  the  y  axis  at  velocity  V2+2c2. 
Also,  on  the  fan  of  straight  characteristics  in  zone  3 

&  =  !L  =  V3-c3  (B.57) 

dt  t 

combining  equations  B.56  and  B.57  to  eliminate  C3  gives 

1H(!+t+c2)  (a88) 

Similarly,  solve  for  c3  by  combining  equations  B.56  and  B.57  to  eliminate  V3  giving 

c3  =  i  (v2  +  2c2  -  I)  (B.59) 

substituting  A3  =  ^C3J  into  equation  B.59  gives 

f)2 

Zone  3  is  bounded  by  characteristics  I  and  II.  For  characteristic  II,  ^  =  V2  -  c2,  and 
for  characteristic  I,  on  the  y  axis,  the  water  boundary  moves  at  V2  +  2 c2  Therefore  the  zone  3 
solution  is  bounded  by 

(V2-c2)t<y  <  (V2  +  2c2)  t  (B.61) 


Calculations  for  Case  V 

Case  V  occurs  when  the  conditions  for  case  III  are  met  but  the  result  for  Ai  in  equation  B.43 
is  less  than  zero.  The  situation  is  illustrated  in  figure  B-8.  To  determine  when  Ai  <  0,  use 
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Figure  B-8:  Case  V,  Fixed  Coordinate  System 


equation  B.43 


which  can  be  rewritten  as 


C2  +  cp  Vi  —  Vo 

rv  ^ 


1  2 


<0 


V<i  +  2c2  <  Vq  —  2cq 


(B.62) 


(B.63) 


The  solution  for  zone  3  in  this  case  is  identical  to  the  solution  for  zone  3  in  case  IV.  The 
solution  for  zone  4  can  be  developed  from  characteristics  similar  to  equations  B.57  through 
B.61.  The  results  are 


U  =  || 

(B.64) 

ii 

3 -*+*»)* 

(B.65) 

The  zone  4  solution  is  bounded  by 


(Vo  -  2co)  t  <  y  <  (Vb  +  cq)  t 


(B.66) 
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Appendix  C 


MATLAB  Program  to  Solve 
Riemann  Problem  Using  the 
Random  Choice  Method 


%  MATLAB  Program  to  Solve  Riemann  Problem  Using 
%  the  Random  Choice  Method 

%  gravity  acceleration  constant 
g=9.8; 

%  wave  form  for  testing 

ymax=.5; 

ymin=-.5; 

dely=.01; 

jmax=(ymax-ymin)/dely; 

%  time  stuff 
delt=0.001; 
tfinal=.08; 
tmax=tfinal/delt ; 

%  Boundary  Conditions  (the  Riemann  Problem  Initial  Conditions) 
vleft=.8; 
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hleft=.6; 

vright=l; 

hright=.2; 

%  set  initial  values 
for  jj=l:jmax 

yplot  (j  j )  =y  min+jj  *  dely-dely  /  2 ; 

h(jj)=hleft; 

v(jj)=vleft; 

if  yplot(jj)>=0 

h(jj)=hright; 

v(jj)=vright; 

end 

end 

for  tt=l:tmax  %  loop  for  time  (t)  evolution 

y=(-l)''tt*dely*rand/2;  %  this  is  the  random  number  for  sampling 
%  each  interval 
%  do  random  choice  1  cycle 

if  y  <=0  %  if  random  number  less  than  zero,  sample  to  left 

for  jj=2:jmax+l 

if  jj  ==  1 

v2=vleft; 

lambda2=hleft; 

vO=v(jj); 

lambdaO==h(jj); 

end 

if  jj  ==  jmax+1  %  Right  hand  side  BC 

v2=v(jj-l); 

lambda2=h(jj-l ) ; 

vO=vright; 

lambdaO=hright ; 
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end 

if  jj  >  1 

if  jj  <  jmax+1 

v2=v(jj-l); 

lambda2=h(jj-l); 

vO=v(jj); 

lambdaO=h(jj); 

end 

end 

%  get  solution  for  the  intervals  Riemann  problem 

soln=riemann(lambda0)lambda2,v0,v2,y,delt,g); 

vnew  (jj- 1 ) =soln  (2) ; 
hnew  (j  j- 1 )  =soln(  1) ; 
end 
end 

if  y  >  0  %  if  random  number  >  0,  sample  to  right 

for  jj=l:jmax 

if  jj  ==  1 

v2=vleft; 

lambda2=hleft ; 

v0=v(jj); 

lambdaO=h(jj); 

end 

if  jj  =  jmax+1 
v2=v(jj-l); 
lambda2=h(jj-l ) ; 
v0=vright; 
lambdaO=hright ; 
end 

if  jj  >  1 
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if  jj  <  jmax+1 

v2=v(jj-l); 

lambda2=h(jj-l); 

vO=v(jj); 

lambdaO=h(jj); 

end 

end 

%  get  solution  for  the  intervals  Riemann  problem 

soln=riemann(lambda0,lambda2,v0,v2,y,delt,g); 

vnew  (jj ) =soln(2) ; 

hnew(jj)=soln(l); 

end 

end 

%  set  up  for  next  time  step 

for  jj=l:jmax 

v(jj)=vnew(jj); 

b(jj)=hnew(jj); 

end 

end  %  end  time  evolution 

%  This  program  returns  height  and  velocity  for  the  Riemann  problem  with 
%  the  initial  conditions  provided  as  calling  arguments.  The  height  and 
%  velocity  are  returned  for  location  y  at  time  t. 
function  soln=riemann(lambda0,lambda2,v0,v2,y,time,g) 

%  soln(l)=height,  soln(2)= velocity 
%  Riemann  Problem  solver 

%  these  are  set  so  that  numbers  are  always  returned 

soln(l)=lambda2; 

soln(2)=v2; 

done=0; 

if  Iambda0==lambda2 
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if  vO  ==  v2 
done=l; 

soln(  1)  =lambdaO; 

soln(2)=v0; 

end 

end 

cO=sqrt(g*lambdaO) ; 
c2=sqrt(g*lambda2) ; 

%  case  4  ******************************************** 

if  lambdaO  <  O.Ol*  lambda2 

done=l; 

soln(l)=(l/(9*g))*(v2+2*c2-y/time)~2; 

soln(2)=(2/3)*(y/time+v2/2+c2); 

if  y  <=  (v2-c2)*time 

soln(  1)  =lambda2; 

soln(2)=v2; 

end 

if  y  >=  (v2+2*c2)*time 

soln(l)=0; 

soln(2)=0; 

end 

end 

*****  en(j  cage  4  ***************************** 

if  done  ==  0 

%  case  3  and  5  *************************************** 
if  (v0-v2)>  2*abs(c0-c2) 
if  (v2+2*c2)>=(vO-2*cO)  %  case  3 
done=l; 

vl=(v0+v2)/2+c2-c0; 

lambdal = ( (v2-v0)  /4  +  (c2+c0)  /  2)  ~  2/g; 
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cl=sqrt(g*lambdal) ; 
if  y  >=(vO+cO)*time 
soln(  1 )  =lambdaO ; 
soln(2)=v0; 
end 

if  y  <=(v2-c2)*time 
soln(  1)  =lambda2; 
soln(2)=v2; 
end 

if  y  >  (v2-c2)*time 

if  y  <=  (vl+cl)*time 
soln(  1)  =lambdal ; 
soln(2)=vl; 
end 

if  y  <=  (vl-cl)*time  %  in  zone  3 

soln(l)=  (1  /  (9*g))*  (v2+2*c2-y/ time)  ~2 

soln(2)=(2/3)*(y/time+v2/2+c2); 

end 

end 

if  y  <  (vO+cO)*time 

if  y  >  (vl+cl)*time  %  in  zone  4 

soln(l)=  (l/(9*g))*(2*c0-v0+y/time)^2 

soln(  2) =(  2/  3)  *  (y/ time+ v2  /  2-cO) ; 

end 

end 

else  %  case  5 
done=l; 

if  y  <=  (v2+2*c2)*time 
soln(l)=(l/ (9*g))*(v2+2*c2-y/time)  "2; 
soln(2) =(2/3)*  (y /timed- v2/2+c2) ; 


if  y  <=  (v2-c2)*time 

soln(  1)  =lambda2 ; 

soln(2)=v2; 

end 

end 

if  y  >=  (v0-2*c0)*time 

soln(l) =(  1  /  (9*g) ) *  (v0-2*c0-y/time) ~  2; 

soln(2) =(2/  3)  *  (y /timed- vO  /  2-cO) ; 

if  y  >=  (vO+cO)*time 

soln(  1)  =lambdaO; 

soln(2)=v0; 

end 

end 

if  y  >  (v2+2*c2)  *time 

if  y  <  (v0-2*c0)*time 

soln(l)=0; 

soln(2)=0; 

end 

end 

end 

end 

%  end  case  3  and  5 
end 

lambdal=0; 
if  done  ==  0 

%  possible  case  1  or  2  ******************************************** 

R=rootfind(vO ,  v2,c0  ,c2) ; 

for  n=l:2 

if  R(n)>0 

r=R(n); 
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root==sqrt  (8*r~  2+1) ; 
zee=cO*r+vO; 

vl=cO*(r-l  /  (4*r)  *  ( 1  +root ) )  +vO; 

cl=c0*sqrt(0.5*(root-l)); 

lambdal=cl  "2^; 

end 

end 

%  test  dillingham  eq  49 
iflambdal  ~=0 

rhs=sqrt(g*  (lambdal  +lambdaO)  *  (lambdal-lambdaO)  ~  2/  (2*lambdal*lambda0) ) ; 
end 

if  (v2-vO)  <  rhs 

vl=c0*(r-l/(4*r)*(l+root))+v0; 

if  y  <=  (v2-c2)*time 

soln(l)=lambda2; 

soln(2)=v2; 

end 

if  y  >  (v2-c2)*time 

if  y  <=  (vl-cl)*time 

soln(  1  ) = ( 1  /  (9*  g) )  *  ( v2+2*c2-y  /  time)  ~  2; 

soln  ( 2) =(2/3)  *  ( (y/time)-vO+ ( v2-v0)  /2+c2) + vO; 

end 

if  y  >  (vl-cl)*time 
if  y  <=  zee*time 
soln(l)  =lambdal ; 
soln(2)=vl; 
end 

if  y  >  zee*time 

soln(l)=lambdaO; 

soln(2)=v0; 
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end 

end 

end 

%  case  2 
else 

%  L  returns  lambda  1  (Dillingham  eq  (55)) 
lambdal=rootfindl  (v0,v2,c0,c2,g); 

vl=vO+sqrt  (g*  (lambdaO  -f  lambda  1 )  *  (lambdal-lambdaO)  ~2/(2*lambda0*lambdal ) ) ; 
zeel=0; 

if  lambdal~=lambdaO 

zeel=(lambdal*vl-lambdaO*vO)  /  (lambdal-lambdaO); 
end 

zee2=0; 

if  lambda2~=lambdal 

zee2=(lambda2*v2dambdal*vl)/ (lambda2-lambdal); 
end 

if  y  >=  zeel*time 
soln(  1)  =lambdaO; 
soln(2)=v0; 
end 

if  y  <  zeel*time 

if  y  >  zee2*time 

soln(  1)  =lambdal ; 

soln(2)=vl; 

end 

end 

if  y  <=  zee2*time 
soln(l)=lambda2; 
soln(2)=v2; 
end 
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end 

%  **  end  case  1  or  2  *********************************************** 
end 
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